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1 Outline: Patterns in the Cosmic Web 

The spatial cosmic matter distribution on scales of a few up to more than a hundred 
Megaparsec displays a salient and pervasive foamlike pattern. Revealed through 
the painstaking efforts of redshift survey campaigns, it has completely revised our 
view of the matter distribution on these cosmological scales. The weblike spatial ar- 
rangement of galaxies and mass into elongated filaments, sheetlike walls and dense 
compact clusters, the existence of large near-empty void regions and the hierarchi- 
cal nature of this mass distribution - marked by substructure over a wide range of 
scales and densities - are three major characteristics of we have come to know as 
the Cosmic Web. 

The vast Megaparsec cosmic web is undoubtedly one of the most striking ex- 
amples of complex geometric patterns found in nature, and the largest in terms 
of sheer size. In a great many physical systems, the spatial organization of mat- 
ter is one of the most readily observable manifestations of the forces and processes 
forming and moulding them. Richly structured morphologies are usually the conse- 
quence of the complex and nonlinear collective action of basic physical processes. 
Their rich morphology is therefore a rich source of information on the combination 
of physical forces at work and the conditions from which the systems evolved. In 
many branches of science the study of geometric patterns has therefore developed 
into a major industry for exploring and uncovering the underlying physics (see e.g. 



Balbus&Hawlev, 1998) 



Computer simulations suggest that the observed cellular patterns are a promi- 
nent a nd natural aspe ct of cosmic structure formation through gravitational insta- 
bility (iPeebles , Il980l) . the standard paradigm for the emergence of structure in our 



Universe. Structure in the Universe is the result of the gravitational growth of tiny 
density perturbations and the accompanying tiny velocity perturbations in the pri- 
mordial Universe. Supported by an impressive body of evidenc e, primarily those of 



temperature fluctuat i ons in the cos r nic rn icrowave background (iSmoot et al.L 11992 



Bennett et al.L 120031 : ISpergel et al.L 120061) . the character of the primordial random 



density and velocity perturbation field is that of a homogeneous and isotropic spatial 
Gaussian process. Such fields of primordial Gaussian perturbations in the gravita- 
tional potential are a natural product of an early inflationary phase of our Universe. 

The early linear phase of pure Gaussian density and velocity perturbations has 
been understood in great depth. This knowledge has been exploited extensively in 
extracting a truely impressive score of global cosmological parameters. Notwith- 
standing these successes, the more advanced phases of cosmic structure formation 
are still in need of substantially better understanding. Mildly nonlinear structures 
do contain a wealth of information on the emergence of cosmic structure at a stage 
features start to emerge as individually recognizable objects. The anisotropic fila- 
mentary and planar structures, the characteristic large underdense void regions and 
the hierarchical clustering of matter marking the weblike spatial geometry of the 
Megaparsec matter distribution are typical manifestations of mildly advanced gravi- 
tational structure formation. The existence of the intriguing foamlike patterns repre- 
sentative for this early nonlinear phase of evolution got revealed by major campaigns 



RA 




Fig. 1. The galaxy distribution uncovered by the 2dF GALAXY REDSHIFT SURVEY. De- 
picted are the positions of 221414 galaxies in the final 2dFGRS catalogue. Clearly visi- 
ble is the foamlike geometry of walls, filaments and massive compact clusters surrounding 
near-empty void regions. Image courtesy of M. Colless; see also Colless et al. (2003) and 
[http ://w w w2 . aao. gov. au/2dFGRS7 1 
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to map the galaxy distribution on Megaparsec scales revealed while ever larger com- 
puter N-body simulations demonstrated that such matter distributions are indeed 
typical manifestations of the gravitational clustering process. Nonetheless, despite 
the enormous progres true insight and physical understanding has remained limited. 
The lack of readily accessible symmetries and the strong nonlocal influences are 
a major impediment towards the development relevant analytical descriptions. The 
hierarchical nature of the gravitational clustering process forms an additional com- 
phcation. While small structures materialize before they merge into large entities, 
each cosmic structure consists of various levels of substructure so that instead of 
readily recognizing one characteristic spatial scale we need to take into account a 
range of scales. Insight into the complex interplay of emerging structures through- 
out the Universe and at a range of spatial scales has been provided through a variety 
of analytical approximations. Computer simulations have provided us with a good 
impression of the complexities of the emerging matter distribution, but for the anal- 
ysis of the resulting patterns and hierarchical substructure the toolbox of descriptive 
measures is still largely heuristic, ad-hoc and often biased in character. 

While cosmological theories are describing the development of structure in 
terms of continuous density and velocity fields, our knowledge stems mainly from 
discrete samplings of these fields. In the observational reality galaxies are the main 
tracers of the cosmic web and it is through measuring the redshift distribution of 
galaxies that we have been able to map its structure. Likewise, simulations of the 
evolving cosmic matter distribution are almost exclusively based upon N-body par- 
ticle computer calculation, involving a discrete representation of the features we 
seek to study. 

Both the galaxy distribution as well as the particles in an N-body simulation are 
examples of spatial point processes in that they are (1) discretely sampled (2) have 
a irregular spatial distribution. The translation of discretely sampled and spatially 
irregularly distributed sampled objects into the related continuous fields is not nec- 
essarily a trivial procedure. The standard procedure is to use a filter to process the 
discrete samples into a representative reconstruction of the underlying continuous 
field. It is the design of the filter which determines the character of the reconstruc- 
tion. 

Astronomical applications are usually based upon a set of user-defined filter 
functions. Nearly without exception the definition of these include pre-conceived 
knowledge about the features one is looking for. A telling example is the use of 
a Gaussian filter. This filter will suppress the presence of any structures on a scale 
smaller than the characteristic filter scale. Moreover, nearly always it is a spherically 
defined filter which tends to smooth out any existing anisotropics. Such procedures 
may be justified in situations in which we are particularly interested in objects of 
that size or in which physical understanding suggests the smoothing scale to be of 
particular significance. On the other hand, they may be crucially inept in situations 
of which we do not know in advance the properties of the matter distribution. The 
gravitational clustering process in the case of hierarchical cosmic structure forma- 
tion scenarios is a particularly notorious case. As it includes structures over a vast 
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Fig. 2. SDSS is the largest and most systematic sky survey in the history of astronomy. It is a 
combination of a sky survey in 5 optical bands of 25% of the celestial (northern) sphere. Each 
image is recorded on CCDs in these 5 bands. On the basis of the images/colours and their 
brightness a million galaxies are subsequently selected for spectroscopic follow-up. The total 
sky area covered by SDSS is 8452 square degrees. Objects will be recorded to ra;„„ = 23.1. In 
total the resulting atlas will contain 10* stars, 10** galaxies and 10^ quasars. Spectra are taken 
of around 10* galaxies, 10' quasars and 10' unusual stars (in our Galaxy). Of the 5 public 
data releases 4 have been accomplished, ie. 6670 square degrees of images is publicly avail- 
able, along with 806,400 spectra. In total, the sky survey is now completely done (107%), the 
spectroscopic survey for 68%. This image is taken from a movie made by Subbarao, Suren- 
dran & Landsberg (see website: http://astro.uchicago.edu/cosmus/projects/sloangalaxies/). It 
depicts the resulting redshift distribution after the 3rd public data release. It concerns 5282 
square degrees and contained 528,640 spectra, of which 374,767 galaxies. 



range of scales and displays a rich palet of geometries and patterns any filter design 
tends to involve a discrimination against one or more - and possibly interesting - 
characteristics of the cosmic matter distribution it would be preferrable to define 
filter and reconstruction procedures that tend to be defined by the discrete point 
process itself. 

A variety of procedures that seek to define and employ more "natural" fil- 
ters have been put forward in recent years. The scale of smoothing kernels can be 
adapted to the particle number density, yielding a density field that retains to a large 
extent the spatial information of the sampled density field. While such procedures 
may still have the disadvantage of a rigid user-defined filter function and filter ge- 
ometry, a more sophisticated, versatile and particularly promising class of functions 
is that of wavelet-defined filters (see contribution B.J. T. Jones). These can be used to 
locate contributions on a particular scale, and even to trace features of a given geom- 
etry. While its successes have been quite remarkable, the success of the application 
is still dependent on the particular class of employed wavelets. 

In this contribution we will describe in extensio the technique of the Delau- 
nay Tessellation Field Estimator. DTFE is based upon the use of the Voronoi 
and Delaunay tessellations of a given spatial point distribution to form the basis 
of a natural, fully self-adaptive filter for discretely sampled fields in which the 
Delaunay tessellations are used as multidimensional interpolation intervals. The 
method has been defined, introduced and developed by ISchaap & van de Wevgaert 
(l2000l) and forms an elaborat ion of the velocity interpolation scheme introduced by 



Bernardeau & van de Wevga ert ( 1996). 



Our focus on DTFE will go along with a concentration on the potential of spatial 
tessellations as a means of estimating and interpolating discrete point samples into 
continuous field reconstructions. In particular we will concentrate on the potential 
of Voronoi and Delaunay tessellations. Both tessellations - each others dual - are 
fundamental concepts in the field of stochastic geometry. Exploiting the character- 
istics of these tessellations, we will show that the DTFE technique is capable of 
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delineating the hierarchical and anisotropic nature of spatial point distributions and 
in outlining the presence and shape of voidlike regions. The spatial structure of the 
cosmic matter distribution is marked by precisely these , and precisely this potential 
has been the incentive for analyzing cosmic large scale structure. DTFE exploits 
three particular properties of Voronoi and Delaunay tessellations. The tessellations 
are very sensitive to the local point density, in that the volume of the tessellation 
cells is a strong function of the local (physical) density. The DTFE method uses this 
fact to define a local estimate of the density. Equally important is their sensitivity 
to the local geometry of the point distribution. This allows them to trace anisotropic 
features such as encountered in the cosmic web. Finally it uses the adaptive and 
minimum triangulation properties of Delaunay tessellations to use them as adap- 
tive spatial interpolation intervals for irregular point distributions. In this it is the 
first order version of the Natural Neighbour method (NN method). The theoretical 
basis for the NN method, a generic smooth and local higher order spatial interpo- 
lation technique developed by ex perts in t he field of c omputational geom etry, has 
been worked out in great detail bv lSibsonl(ll980i 1 198 lb and WatsonI (119921). As has 
been demonstrated b y telling example s in geophysics (IBraun & Sambridgei 119951) 
and solid mechanics (iSukuman Il998h NN methods hold tremendous potential for 
grid-independent analysis and computations. 

Following the definition of the DTFE technique, we will present a systematic 
treatment of various virtues of relevance to the cosmic matter distribution. The re- 
lated performance of DTFE will be illustrated by means of its success in analyzing 
computer simulations of cosmic structure formation as well as that of the galaxy 
distribution in large-scale redshift surveys such as the 2dFGRS and SDSS surveys. 
Following the definition and properties of DTFE, we will pay attention to extensions 
of the project. Higher order Natural Neighbour renditions of density and velocity 
fields involve improvements in terms of smoothness and discarding artefacts of the 
underlying tessellations. 

Following the determination of the "raw" DTFE produces density and veloc- 
ity fields - and/or other sampled fields, e.g. temperature values in SPH simula- 
tions - the true potential of DTFE is realized in the subsequent stage in which 
the resulting DTFE fields get processed and analyzed. We will discuss an ar- 
ray of techniques which have been developed to extract information on aspects 
of the Megaparsec matter distribution. Straightforward processing involves sim- 
ple filtering and the production of images from the reconstructed field. Also, 
we will shortly discuss the use of the tessellation fields towards defining new 
measures for the topology and geometry of the underlying matter distribution. 
The determination of the Minkowski functionals of isodensity surfac e s, following 



the SURFGEN formahsm dSahni. Sathvaprakash & ShandarinI Il998l: Isheth et al- 
20031: IShandarin. Sheth & SahniL |2004 . can be greatly facilitated on the basis of 



the Delaunay triangulation itself and the estimated DTFE density field. Extending 
the topological information contained in the Minkowski functionals leads us to the 
concept of Betti numbers and a-shapes, a filtration of the Delaunay complex of 



a dataset (lEdelsbrunner. Kirkpatrick & Seidell 1 1 9831 : lEdelsbrunner & Muckelll994l: 
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Edelsbrunner. Letscher & ZomorodianI l2002h . With the associated persistence di- 



agrams the a-shapes encode the evolution of the Betti numbers. We have started 
to add to the arsenal of t ools t o quantify the pa tterns in the Megaparsec matter 
distribution (iVegter et al. , I2OO7I: iEldering, l2006h . Particularly interesting are the 
recentl y developed elaborate and advanced tec hniques of watershed void identifi- 
cation (jPlaten. van de W evgaert & JonesL 2007 ) and Multiscale Morphology Filter 
( Aragon-Calvol 20071 : Aragon-C alvo et all " 2007 ). These methods enable the unbi- 
ased identification and measurement of voids, walls, filaments and clusters in the 
galaxy distribution. 

Preceding the discussion on the DTFE and related tessellation techniques, we 
will first have to describe the structure, dynamics and formation of the cosmic web. 
The complex structure of the cosmic web, and the potential information it does 
contain, has been the ultimate reason behind the development of the DTFE. 



2 Introduction: The Cosmic Web 

Macroscopic patterns in nature are often due the collective action of basic, often 
even simple, physical processes. These may yield a surprising array of complex 
and genuinely unique physical manifestations. The macroscopic organization into 
complex spatial patterns is one of the most striking. The rich morphology of such 
systems and patterns represents a major source of information on the underlying 
physics. This has made them the subject of a major and promising area of inquiry. 



2.1 Galaxies and the Cosmic Web 



One of the most striking examples of a physical system displaying a salient geo- 
metrical morphology, and the largest in terms of sheer size, is the Universe as a 
whole. The past few decades have revealed that on scales of a few up to more than a 
hundred Megaparsec, the galaxies conglomerate into intriguing cellular or foamlike 
patterns that pervade throughout the observable cosmos. An initial hint of this cos- 
mic web was seen in the view of the local U niverse offered by the first Cf A redshift 
slice (Ide Lapparent. Geller & Huchralll986L . In recent years this view has been ex- 
panded dramatically to the present grand vistas offered by the 100,00 0s of galaxies 
in the 2dF - two-degr ee field - Galaxy Reds hift Survey, the 2dFGRS ( Colless et al 



2003 ) and SDSS (e.g. Tegmark et al. , 2004) galaxy redshift surveys.!^. Galaxies are 
found in dense, compact clusters, in less dense filaments, and in sheetlike walls 
which surround vast, almost empty regions called voids. This is most dramatically 
illustrated by the 2dFGRS and SDSS maps. The published maps of the distribution 
of nearly 250,000 galaxies in two narrow "slice" regions on the sky yielded by the 
2dFGRS surveys reveal a far from homogeneous distribution (fig. [1). Instead, we 



^See |http : //www . mso . anu . edu . au/2dFGRS/ | and |http : //www . sdss . org/ 1 



Geometric Analysis Cosmic Web 9 




Fig. 3. Equatorial view of the 2MASS galaxy catalog (6h RA at centre). The grey-scale rep- 
resents the total integrated flux along the line of sight - the nearest (and therefore brightest) 
galaxies produce a vivid contrast between the Local Supercluster (centre-left) and the more 
distant cosmic web. The dark band of the Milky Way clearly demonstrates where the galaxy 
catalog becomes incomplete due to source confusion. Some well known large-scale structures 
are indicated: P-P=Perseus-Pisces supercluster; H-R=Horologium-Reticulum supercluster; 
P-I=Pavo-Indus supercluster; GA= 'Great Attractor'; GC=Galactic Centre; S-C=Shapley 
Concentration; 0-C=Ophiuchus Cluster; Virgo, Coma, and Hercules=Virgo,Coma and Her- 
cules superclusters. The Galactic 'anti-centre' is front and centre, with the Orion and Taurus 
Giant Molecular Clouds forming the dark circular band near the centre. Image courtesy of 
J.H. Jarrett. Reproduced with permission from the Publications of the Astronomical Society 
of Australia 21(4): 396-403 (T.H. Jarrett). Copyright Astronomical Society of Australia 2004. 
Published by CSIRO PUBLISHING, Melbourne Australia. 



recognize a sponge-like arrangement, with galaxies aggregating in striking geomet- 
ric patterns such as prominent filaments, vaguely detectable walls and dense com- 
pact clusters on the periphery of giant voido The three-dimensional view emerging 
from the SDSS redshift survey provides an even more convincing image of the intri- 
cate patterns defined by the cosmic web (fig.[I]l. A careful assessment of the galaxy 
distribution in our immediate vicintiy reveals us how we ourselves are embedded 
and surrounded by beautifully delineated and surprisingly sharply defined weblike 

^It is important to realize that the interpretation of the Megaparsec galaxy distribution is 
based upon the tacit yet common assumption that it forms a a fair reflection of the underly- 
ing matter distribution. While there are various indications that this is indeed a reasonable 
approximation, as long as the intricate and complex process of the formation of galaxies has 
not been properly understood this should be considered as a plausible yet heuristic working 
hypothesis. 
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Fig. 4. The CfA Great Wall (bottom slice, Geller & Huchra 1989) compared with the Sloan 
Great Wall (top slice). Both structures represent the largest coherent structural in the galaxy 
redshift surveys in which they were detected, the CfA redshift survey and the SDSS redshift 
survey. The (CfA) Great Wall is a huge planar concentration of galaxies with dimensions 
that are estimated to be of the order of 6O/2"' x 170/!"' x 5/z"' Mpc. Truely mindboggling 
is the Sloan Great Wall, a huge conglomerate of clusters and galaxies. With a size in the 
order of 400/?" 'Mpc it is at least three times larger than the CfA Great Wall. It remains to 
be seen whether it is a genuine physical structure or mainly a stochastic arrangement and 
enhancement, at a distance coinciding with the survey's maximum in the radial selection 
function. Image courtesy of M. Juric, see also Gott et al. 2005. Reproduced by permission of 
the AAS. 



Structures. In particular the all-sky nearby infrared 2MASS survey ( see fig. |3]l pro - 
vides us with a meticulously clear view of the web surrounding us d Jarre tti 120041) . 
The cosmic web is outlined by galaxies populating huge filamentary and wall-like 
structures, the sizes of the most conspicuous one frequently exceeding lOOh ' Mpc. 
The closest and best studied of these massive anisotropic matter concentrations can 
be identified with known supercluster complexes, enormous structures comprising 
one or more rich clusters of galaxies and a plethora of more modestly sized clumps 
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Fig. 5. The cosmic web at high redshifts: a prominent weblike features at a redshift z ~ 
3. 1 found in a deep view obtained by the Subaru telescope. Large scale sky distribution of 
283 strong Lyo- emitters (black filled circles), the Lytt absorbers (red filled circles) and the 
extended hya emitters (blue open squares). The dashed lines indicate the high-density region 
of the strong Lyo- emitters. From Hayashino et al. 2004. Reproduced by permission of the 
A AS. 
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of galaxies. A prominent and representative nearby specimen is the Perseus-Pisces 
supercluster, a 5/;"' wide ridge of at least Mpc length, possibly extending out 
to a total length of 140/?"' Mpc. While such giant elongated structures are amongst 
the most conspicuous features of the Megaparsec matter distribution, filamentary 
features are encountered over a range of scales and seem to represent a ubiquitous 
and universal state of concentration of matter In addition to the presence of such fil- 
aments the galaxy distribution also contains vast planar assemblies. A striking local 
example is the Great Wall, a huge planar concentration of galaxies with dimensions 
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that a re estimated to be of the order of 60h ' x 1 70h ' x 5/; ' Mpc JOeller & Huchra , 
19891) . In both the SDSS and 2dF surveys even more impressive planar complexes 
were recognized, with dimensions substantially in excess of those of the local Great 
Wall. At the moment, the socalled SDSS Great Wall appears to be the largest known 
structure in the Universe (see fig.|4]i. Gradually galaxy surveys are opening the view 
onto the large scale distribution of galaxies at high redshifts. The Subaru survey 
has even managed to map out a huge filamentary feature at a redshift of z ~ 3.1, 
perhaps the strongest evidence for the existence of pronounced cosmic structure at 
early cosmic epochs (fig.|5]l. 




Fig. 6. Comparison of optical and X-ray images of Coma cluster. Top: optical image. Image 
courtesy of O. Lopez-Cruz. Lower: X-ray image (ROSAT). 



Of utmost significance for our inquiry into the issue of cosmic structure forma- 
tion is the fact that the prominent structural components of the galaxy distribution 
- clusters, filaments, walls and voids - are not merely randomly and independently 
scattered features. On the contrary, they have arranged themselves in a seemingly 
highly organized and structured fashion, the cosmic foam or cosmic web. They are 
woven into an intriguing/oam/ife tapestry that permeates the whole of the explored 
Universe. The vast under-populated void regions in the galaxy distribution represent 
both contrasting as well as complementary spatial components to the surrounding 
planar filamentary density enhancements. At the intersections of the latter we 
often find the most prominent density enhancements in our universe, the clusters of 
galaxies. 
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2.2 Cosmic Nodes: Clusters 



Within and around these anisotropic features we find a variety of density condensa- 
tions, ranging from modest groups of a few galaxies up to massive compact galaxy 
clusters. The latter stand out as the most massive, and most recently, fully collapsed 
and virialized objects in the Universe. Approximately 4% of the mass in the Uni- 
verse is assembled in rich clusters. They may be regarded as a particular population 
of cosmic structure beacons as they typically concentrate near the interstices of the 
cosmic web, «ot/e^ formi ng a recognizable tracer of the cosmic matter distribution 
dBorgani & GuzzoluOOlh . Clusters not only function as wonderful tracers of struc- 
ture over scales of dozens up to hundred of Megaparsec (fig. |2.2t but also as useful 
probes for precision cosmology on the basis of their unique physical properties. 

The richest clusters contain many thousands of galaxies within a relatively small 
volume of only a few Megaparsec size. For instance, in the nearby Virgo and Coma 
clusters more than a thousand galaxies have been identified within a radius of a 
mere 1.5/i ' Mpc around their core (see fig. 12.1b . Clusters are first and foremost 
dense concentrations of dark matter, representing overdensities A ~ 1000. In a sense 
galaxies and stars only form a minor constituent of clusters. The cluster galaxies are 
trapped and embedded in the deep gravitational wells of the dark matter. These are 
identified as a major source of X-ray emission, emerging from the diffuse extremely 
hot gas trapped in them. While it fell into the potential well, the gas got shock-heated 
to temperatures in excess of T > 10^ K, which results in intense X-ray emission due 
to the bremsstrahlung radiated by the electrons in the highly ionized intracluster 
gas. In a sense clusters may be seen as hot balls of X-ray radiating gas. The amount 
of intracluster gas in the cluster is comparable t o that locked into stars, and stands 
for QjcM ~ 0.0018 (iFukugita & Peeblesl 120041) . The X-ray emission represents a 
particularly useful signature, an objective and clean measu re of the potential well 
depth , directly related to the total mass of the cluster (see e.g. lReiprich & Bohringer , 
1999h . Through their X-ray brightness they can be seen out to large cosmic depths. 
The deep gravitational dark matter wells also strongly affects the path of passing 
photons. While the resulting strong lensing arcs form a spectacular manifestation, 
it has been the more moderate distortion of backgroun d galaxy images in the weak 
lensing regime (iKaiserl Il992t iKaiser & Squiresl Il993h which has opened up a new 
window onto the Universe. The latter has provided a direct probe of the dark matter 



content of cluste rs and the large scale universe (for a review see e.g. lMelUeiill999 
Refregieiil2003l) . 



2.3 Cosmic Depressions: the Voids 

Complementing this cosmic inventory leads to the existence of large voids, enor- 
mous regions with sizes in the range of 20 - SO/z"' Mpc that are practically devoid 
of any galaxy, usually roundish in shape and occupying the major share of space 
in the Universe. Forming an essential ingredient of the Cosmic Web, they are sur- 
rounded by elongated filaments, sheetlike walls and dense compact clusters. 
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Fig. 7. The spatial cluster distribution. The full volume of the X-ray REFLEX cluster survey 
within a distance of 600h"'Mpc. The REFLEX galaxy cluster catalogue (Bohringer et al. 
2001), contains all clusters brighter than an X-ray flux of 3x 10^'^ergs^'cm"- over a large part 
of the southern sky. The missing part of the hemisphere delineates the region highly obscured 
by the Galaxy. Image courtesy of S. Borgani & L. Guzzo, see also Borgani & Guzzo (2001). 
Reproduced by permission of Nature. 



Voids have been known as a feature of galaxy surveys since the first surveys were 



comp iled ( IChincarini & Roodlll97 5':'Gregorv & Thompson'.'igTS'/ Einasto. Joeveer & Saar 



19801) . Following the discovery bv lKirshner et al.. C198L.1987.) of the most outstand 
ing specimen, the Bootes void, a hint of their cent ral position within a weblike ar- 
range ment came with the first CfA redshift slice dde Lapparent. Geller & Huchra , 
198d) . This view has been expanded dramatically as ma ps of the spatial dist ribu- 



tion of hundreds of tho usands of galaxi e s in th e 2dFGRS jColless et all l2003h and 



SDSS redshift survey jAbazaiian et al.L 120031) became available, recently supple 



mented with a high-resolution study of voids in the nearby Universe based upon 



Geometric Analysis Cosmic Web 



15 




Fig. 8. A region of the 6dF redshift survey marked by the presence of various major voids. 
The image concerns a 3D rendering of the galaxy distribution in a 1000 km/s thick slice along 
the supergalactic SGX direction, at SGX=-2500 km/s. Image courtesy of A. Fairall. 
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the 6dF survey (iHeathJonesl 120041: iFairalM, l2007h . The 2dFGRS and SDSS maps of 
figs. [Hand [1] and the void map of the 6dF survey in fig. |2.2| form telHng illustrations. 

V oids in the galaxy distribution accoun t for about 95% of the total volume (see 
Kauffmann & Fai rall, 19 91; El-Ad, Piran & da Costal Il996l: IlkAd & Piranl Il997 ; 



Hovle & Vogelevl I2OO2I: IPlionis & Basilakosll2002l:lRoias et al.Ll2005b . The tvpical 



sizes of voids in the galaxy distribution depend on the galaxy population used to 
define the voids. Voids defined by galaxies brighter than a typical L* galaxy tend 
to have diameters of order 10 - 20/i"'Mpc, but voids associated with rare luminous 
galaxies can be consid erably larger; diam eters in th e range of 20/;'' - 50/ i~'Mpc 
are not uncommon (e.g. lHovle & Vogele 3 [2002; Pli onis & Basilakosll2002l) . These 
large sizes mean that only now we are beginning to probe a sufficiently large cos- 
mological volume to allow meaningful statistics with voids to be done. As a result, 
the observations are presently ahead of the theory. 



3 Cosmic Structure Formation: 

from Primordial Noise to Cosmic Web 

The fundamental cosmological importance of the Cosmic Web is that it comprises 
features on a typical scale of tens of Megaparsec, scales at which the Universe still 
resides in a state of moderate dynamical evolution. Structures have only freshly 
emerged from the almost homogeneous pristine Universe and have not yet evolved 
beyond recognition. Therefore they still retain a direct link to the matter distribution 
in the primordial Universe and thus still contain a wealth of direct information on 
the cosmic structure formation process. 

In our exploration of the cosmic web and the development of appropriate tools 
towards the analysis of its structure, morphology and dynamics we start from the the 
assumption that the cosmic web is traced by a population of discrete objects, either 
galaxies in the real observational world or particles in that of computer simulations. 
The key issue will be to reconstruct the underlying continuous density and velocity 
field, retaining the geometry and morphology of the weblike structures in all its 
detail. 

In this we will pursue the view that filaments are the basic elements of the cosmic 
web, the key features around which most matter will gradually assemble and the 
channels along which matter is transported towards the highest density knots within 
the network, the clusters of galaxies. Likewise we will emphasize the crucial role of 
the voids - the large underdense and expanding regions occupying most of space - 
in the spatial organization of the various structural elements in the cosmic web. One 
might even argue that it are the voids which should be seen as the key ingredients of 
the cosmic matter distribution. This tends to form the basis for geometrical models 
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of the Megaparsec sc ale matter distribution, with the Voro noi model as its main 
representative (see e.g lvan de Wevgaertlll99lll2002ll2007allbP l 



3.1 Gravitational Instability 

The standard pa radigm of cosmic structure form ation is that of gravititional insta- 
bility scenarios ( Peebles , 1980l : Zeldovich , 1970l) . Structure in the Universe is the 
result of the gravitational growth of tiny density perturbations and the accompany- 
ing tiny velocity perturbations in the primordial Universe. Supported by an impres- 
sive body of evidenc e, primarily that of t emperature fluctuations in the cosmic mi - 



crowave background (ISmoot et al.L Il992t iBennett et al.L l2003t ISpergel et al.L l2006h 



the character of the primordial random density and velocity perturbation field is that 
of a homogeneous and isotropic spatial Gaussian process. Such fields of primordial 
Gaussian perturbations in the gravitational potential are a natural product of an early 
inflationary phase of our Universe. 

The formation and moulding of structure is a result to the gravitational growth 
of the primordial density- and velocity perturbations. Gravity in slightly overdense 
regions will be somewhat stronger than the global average gravitational decelera- 
tion, as will be the influence they exert over their immediate surroundings. In these 
regions the slow-down of the initial cosmic expansion is correspondingly stronger 
and, when the region is sufficiently overdense it may even come to a halt, turn around 
and start to contract. If or as long as pressure forces are not sufficient to counteract 
the infall, the overdensity will grow without bound, assemble more and more matter 
by accretion of matter from its surroundings, and ultimately fully collapse to form 
a gravitationally bound and virialized object. In this way the primordial overdensity 
finally emerges as an individual recognizable denizen of our Universe, their precise 
nature (galaxy, cluster, etc.) and physical conditions determined by the scale, mass 
and surroundings of the initial fluctuation. 



3.2 Nonlinear Clustering 

Once the gravitational clustering process has progressed beyond the initial linear 
growth phase we see the emergence of complex patterns and structures in the density 
field. Highly illustrative of the intricacies of the structure formation process is that 
of the state-of- t he-art N-body computer simulation, the Millennium simulation by 



dSpringel et al 



2005b . Figure |9] shows four time frames out of this massive 10'" 
particle simulation of a AC DM matter distribution in a SOO/i 'Mpc box. The time 
frames correspond to redshifts z - 8.55, z - 5.72, z = 1.39 and z = (ie. at epochs 
600 Myr, 1 Gyr, 4.7 Gyr and 13.6 Gyr after the Big Bang). The earlies time frame is 
close to that of the condensation of the first stars and galaxies at the end of the Dark 
Ages and the reionization of the gaseous IGM by their radiation. The frames contain 



'These Voronoi models are spatial models for cellular/weblike galaxy distributions, not to 
be confused with the application of Voronoi tessellations in DTFE and the tessellation based 
methods towards spatial interpolation and reconstruction 
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Fig. 9. The Cosmic Web in a box: a set of four time slices from the Millennium simulation 
of the ACDM model. The frames show the projected (dark) matter distribution in slices of 
thickness 15/r'Mpc, extracted at z = 8.55, z = 5.12, z = 1.39 and z = 0. These redshifts 
correspond to cosmic times of 600 Myr, 1 Gyr, 4.7 Gyr and 13.6 Gyr after the Big Bang. 
The four frames have a size of 125/7"'Mpc. The evolving mass distribution reveals the major 
characteristics of gravitational clustering: the formation of an intricate filamentary web, the 
hierarchical buildup of ever more massive mass concentrations and the evacuation of large 
underdense voids. Image courtesy of V. Springel & Virgo consortium, also see Springel et al. 
2005. 

the Dark Matter particle distribution in a IS/i 'Mpc thick slice of a 125/i 'Mpc 
region centered on the central massive cluster of the simulation. 

The four frames provide a beautiful picture of the unfolding Cosmic Web, start- 
ing from a field of mildly undulating density fluctations towards that of a pro- 
nounced and intricate filigree of filamentary features, dented by dense compact 
clumps at the nodes of the network. Clearly visible is the hierarchical nature in 
which the filamentary network builds up. At first consisting of a multitude of small 
scale edges, they quickly merge into a few massive elongated channels. 

Large N-body simulations like the Millennium simulation and the many others 
currently available all reveal a few "universal" characteristics of the (mildly) nonUn- 
ear cosmic matter distribution. Three key characteristics of the Megaparsec universe 
stand out: 

• Hierarchical clustering 

• Weblike spatial geometry 

• Voids 

Hierarchical clustering implies that the first objects to condense first are small and 
that ever larger structures form through the gradual merging of smaller structures. 
Usually an object forms through the accretion of all matter and the fusion of all sub- 
structure within its realm, including that of the small-scale objects which had con- 
densed out at an earlier stage. The second fundamental aspect is that of anisotropic 
gravitational collapse. Aspherical overdensities, on any scale and in any scenario, 
will contract such that they become increasingly anisotropic. At first they turn into 
a flattened pancake, rapidly followed by contraction into an elongated filament and 
possibly, dependent on scale, total collapse into a galaxy or a cluster may follow. 
This tendency to collapse anisotropically finds its origin in the intrinsic primordial 
flattening of the overdensity, augmented by the anisotropy of the gravitational force 
field induced by the external matter distribution (i.e. by tidal forces). It is evidently 
the major agent in shaping the weblike cosmic geometry. The third manifest feature 
of the Megaparsec Universe is the marked and dominant presence of large roundish 
underdense regions, the voids. They form in and around density troughs in the pri- 



20 



Rien van de Weygaert & Willem Schaap 



mordial density field. Because of their lower interior gravity they will expand faster 
than the rest of the Universe, while their internal matter density rapidly decreases as 
matter evacuates their interior. They evolve in the nearly empty void regions with 
sharply defined boundaries marked by filaments and walls. Their essential role in 
the organization of the cosmic matter distribution got recognized early after their 
discovery. Recently, their emergence and ev olution has been e xplained within the 
context of hierarchical gravitational scenarios ISheth & van de W eygaert (l2004h . 

The challenge for any viable analysis tool is to trace, highlight and measure 
each of these aspects of the cosmic web. Ideally it should be able to do so without 
resorting to user-defined parameters or functions, and without aff'ecting any of the 
other essential characteristics. We will argue in this contributation that the DTFE 
method, a linear version of natural neighbour interpolation, is indeed able to deal 
with all three aspects (see fig.fTTTi. 



4 Spatial Structure and Pattern Analysis 

Many attempts to describe, let alone identify, the features and components of the 
Cosmic Web have been of a mainly heuristic nature. There is a variety of statistical 
measures characterizin g specific aspects of the large scale matter distribution (for 
an extensive review see lMartmez & Saan,l2002h . For completeness and comparison. 



we list briefly a selection of methods for structure characterisation and finding. It is 
perhaps interesting to note two things about this list: 

a) each of the methods tends to be specific to one particular structural entity 

b) there are no explicit wall-finders. 

This emphasises an important aspect of our Scale Space approach: it provides a 
uniform approach to finding Blobs, Filaments and Walls as individual objects that 
can be catalogued and studied. 



4.1 Structure from higher moments 

The clustering of galaxies and matter is most c ommonly described in terms of a hier- 
archy of correlation functions dPeebles , Il980l) . The two-p oint correlation function - 
and the equiva l ent po wer spectrum, its Fourier transform jPeacock & Doddsl 1 1994 



Tegmark et al.L l2004h - remains the mainstay of cosmological clustering analysis 



and has a solid physical basis. However, the nontrivial and nonlinear patterns of the 
cosmic web are mostly a re s ult of t he phase correlat i ons in the cosmic matter distri - 
bution dRyden & GramannL Il99ll: IChiang & Colesi I2OO0I: IColes & Chia ng. 200i 
While this information is contained in the moments of c ell counts (IPeebles. .198 



de Lapparent. Geller & Huchrai I199II; iGaztanagai 119921) and, more formally so, in 



1X1- 



the full hierarchy of M-point correlation functions ^m, exc ept for the lowest orders 
their measurement has proven to be practically unfeasible (iPeebles , Il98(]l: l ISzapudiL 



Geometric Analysis Cosmic Web 



21 



1998tlJones et al.U2005h . Problem remains that these higher order correlation func- 



tions do not readily translate into a characterization of identifiable features in the 
cosmic web. 



The Void probability Function ( IWhitelll979l : lLachieze-Rev. da Costa & Maurogordato , 



19921) provides a characterisation of the "voidness" of the Universe in terms of a 



function that combined information from many higher moments of the point distri- 
bution. But, again, this did not provide any identification of individual voids. 



4.2 Topological methods 



The shape of the local matter distribution may be trace d on the basis of an anal- 
ysis of the statis tical proper ties of its inertial momen ts ("Bab ul &Starkmani 1 19921: 
Luo & Vishniad [7995; Basi lakos. Plionis & Rowan-R obinson, 2001). These con- 
cepts are closely related to the full characterizatio n of the topology of the matter dis- 
tribution in terms of four Minkowski functionals (iMecke. Buchert & Wagnen.ll994c 
Schmalzing et al 1E999). They are solidly based on the theory of spatial statistics 
and also have the great advantage of being known analytically in the case of Gaus- 
sian random fields. In particular, the genus of the density field has received sub- 
stantial attentio n as a strongly discriminat i ng factor between intrinsically different 
spatial patterns dOott. Dickinson & Melotill986HHoyle & Vogeleyll2002l) . 

The Minkowski functionals provide global characterisations of structure. An 
attempt to extend its scope towards providing locally defined topological mea- 
sures of the density field has been develope d in the SURFGEN project defined 
by Sahni and Shandarin and their coworkers dSahni. Sathvaprakash & Shandarin , 



1998t ISheth et all 120031: IShandarin. Sheth & SahniL 12004 . The main problem re 



mains the user-defined, and thus potentially biased, nature of the continuous density 
field inferred from the sample of discrete objects. The usual filtering techniques 
suppress substructure on a scale smaller than the filter radius, introduce artificial 
topological features in sparsely sampled regions and diminish the flattened or elon- 
gated morphology of the spatial patterns. Quite possibly the introduction of more 
advanced geometry based methods to trace the den s ity fie ld may prove a major ad- 
vance towards solving this problem. iMartmez et al.l (120051) have generalized the use 
of Minkowski Functionals by calculating their values in variously smoothed volume 
limited subsamples of the 2dF catalogue. 



4.3 Cluster and Filament finding 

In the context of analyzing distributions of galaxies we can think of cluster finding 
algorithms. There we might define a cluster as an aggregate of neighbouring galax- 
ies sharing some localised part of velocity space. Algorithms like HOP attempt to 
do this. However, there are always issues arising such as how to deal with substruc- 
ture: that perhaps comes down to the definition of what a cluster is. Nearly always 
coherent structures are identified on the basis of particle positions alone. Velocity 
space data is often not used since there is no prior prejudice as to what the velocity 
space should look like. 
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The connectedness of elongated supercluster structures in the cosmic matter dis- 
tribution was first probed by means of percolation analysis, introduced and empha- 
sized by Zel'dovich and coworkers dZeldovich. Einasto & Shandarini Il982h . while 
a related graph-theoretical construct, the minimum spanning tree of the galaxy 
distribution, was extensively probed and analyse d by Bhavsar and collaborators 
( Barrow. Bhavsar & Sonodal 1985 : Colberg , 2007 ) in an attempt to develop an ob- 
jective measure of filamentarity. 

Finding filaments joi ning neighbouring clusters has been t ackled, u sing quite 
differe nt techniques, by IColberg. Krughoff & ConnoUvl (l2005b and by iPimbblet 
(120051) . More general filament finders h ave been put forward by a number of authors. 
Skeleton analysis of the density field dNovikov. Colombi & Dora, 120061) describes 
continuous density fields by relating density field gradients to density maxima and 
saddle points. This is computationally intensive but quite effective, though it does 
de pend on the a r tefact s in the reconstruction of the continuous density field. 



Stoica et al.l (120051) use a generalization of the classical Candy model to locate 



and catalogue filaments in galaxy surveys. This appraoch has the advantage that it 
works directly with the original point process and does not require the creation of a 

continuous density field. However, it is very computationally intensive^ 

A recently developed method , the Multiscale Morphology Filter ( Aragon-Calvoi 



2007t lAragon-Calvo et al.l l2007h (see sect. 115.6b . seeks to identify different mor- 



phological features over a range of scales. Scale space analysis looks for structures 
of a mathematically specified type in a hierarchical, scale independent, manner. It 
is presumed that the specific structural characteristic is quantified by some appro- 
priate parameter (e.g.: density, eccentricity, direction, curvature components). The 
data is filtered to produce a hierarchy of maps having different resolutions, and at 
each point, the dominant parameter value is selected from the hierarchy to construct 
the scale independent map. While this sounds relatively straightforward, in practice 
a number of things are required to execute the process. There must be an unam- 
biguous definition of the s tructure-defining characteristic. The implementation of 



Aragon-Calvo et al.l (120071) uses the principal components of the local curvature of 



the density field at each point as a morphology type indicator. This requires that 
the density be defined at all points of a grid, and so there must be a method for 
going from a discrete point set to a grid sampled continuous density field. This is 
done using the DTFE methodology since that does minimal damage to the structural 
morphology of the density field (see sect. |8]l. 



4.4 Void Finding 

Voids are distinctive and striking features of the cosmic web, yet finding them sys- 
tematically in surveys and simulations has proved rather difficult. There have been 
extensive searches for voids in galaxy catalogues and in numerical simulations (see 
sect. 12.3b . Identifying voids and tracing their outline within the complex spatial ge- 
ometry of the Cosmic Web appears to be a nontrivial issue. The fact that voids are 
almost empty of galaxies means that the sampl i ng den sity plays a key role in deter- 
mining what is or is not a void (ISchmidt et all 120011) . There is not an unequivocal 
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definition of what a void is and as a resu lt there is considerable disagreement on the 



precise outline of such a region (see e.g. lShandarin et al.Ll2006h . 



Moreov er, void finders are often pre dicated on building vo id structures out o f 



cubic cells ( iKaufiinann & Fairal]lll991l) or out of spheres (e.g. lPatiri et al.L 120061) . 



Such methods attempt to synthesize voids from the intersection of cubic or spheri- 
cal elements and do so with varying degrees of success. Because of the vague and 
different definitions, and the ran ge of different interests in voids, there is a plethora 
of void identification procedures (IKauffmann & Fairalll |199UIE1-Ad & PiranLll997l: 
Aikio & Mah6nenl'l998';'Hovle & Vogelev', '20021; lArbabi-Bidgoli & Miilleil 



2002; 



Plionis & Basilakos, 2002; Patiri et al, 2006; Colberg et al., 2005b; Shandar in et al 



20061) . The "voidfinder" algorithm of isi-Ad & Piranl7l997l) has been at the basis of 
most voidfinding methods. However, this succesfull approach will not be able to an- 
alyze complex spatial configurations in which voids may have arbitrary shapes and 
co ntain a range and variety of substructures. The Void Finder Comparison Project 
of lColberg. Pearce et al.l ( 120071) will clarify ma ny of these issues. 



The watershed-based WVF algorithm of Pl aten, van de Wevgaert & Jonesl(l2007 



see sect. 115.71 1 aims to avoid issues of both sampling density and shape. This new 
and objective voidfinding formalism has been specifically designed to dissect in a 
selfconsistent manner the multiscale character of the void network and the web- 
like features marking its bo undaries. The Watershed Void Finder (WVF) is based 
on the watershed algorith m (iBeucher & Lantuejouli 119791; iMever & Beucher , 1 19901: 
Beucher & Meyen.ll993h . This is a concept from the field of mathematical morphol- 
ogy and image analysis. The WVF is defined with respect to the D TFE density field 
of a discrete point distribution (iSchaap & van de Weygaertll2000l) . assuring optimal 
sensitivity to the morphology of spatial structures and an unbiased probe over the 
full range of substructure in the mass distribution. Because the WVF void finder 
does not impose a priori constraints on the size, morphology and shape of a voids it 
has the potential to analyze the intricacies of an evolving void hierarchy. 



5 Structural Reconstruction 

In the real world it is impossible to get exhaustive values of data at every desired 
point of space. Also astronomical observations, physical experiments and computer 
simulations often produce discretely sampled data sets in two, three or more di- 
mensions. This may involve the value of some physical quantity measured at an 
uregularly distributed set of reference points. Also cosmological theories describe 
the development of structure in terms of continuous (dark matter) density and ve- 
locity fields while to a large extent our knowledge stems from a discrete sampling 
of these fields. 

In the observational reality galaxies are the main tracers of the cosmic web and 
it is mainly through measuring of the redshift distribution of galaxies that we have 
been able to map its structure. Another example is that of the related study of cos- 
mic flows in the nearby Universe, based upon the measured peculiar velocities of 
a sample of galaxies located within this cosmic volume. Likewise, simulations of 
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the evolving cosmic matter distribution are almost exclusively based upon N-body 
particle computer calculation, involving a discrete representation of the features we 
seek to study. Both the galaxy distribution as well as the particles in an N-body 
simulation are examples of spatial point processes in that they are 

discretely sampled 

have an irregular spatial distribution. 

The principal task for any formalism seeking to process the discretely sampled 
field is to optimally retain or extract the required information. Dependent on the 
purpose of a study, various different strategies may be followed. One strategy is 
to distill various statistical measures, or other sufficiently descriptive cosmological 
measures, characterizin g specific aspects of the large scale matter distribution (see 
Martinez & Saan, l2002l also see sect. |4|i. In essence this involves the compression 
of the available information into a restricted set of parameters or functions, with 
the intention to compare or relate these to theoretical predictions. The alternative 
is to translate the discretely sampled and spatially irregularly distributed sampled 
objects into related continuous fields. While demanding in itself, it is complicated 
by the highly inhomogeneous nature of the sample point distribution. The transla- 
tion is a far from trivial procedure. If forms the subject of an extensive literature 
in computer science, visualization and applied sciences. An interesting comparison 
and application of a few different techniques is shown in fig. [El It shows how some 
methods to be discussed in the following sections fare when applied to the recon- 
struction of Martian lands capes from measurements by the MOLA instrument on 
the Mars Global Surveyor jAbramov & McEwenll2004 . 



5.1 Spatial Data: Filtering and Interpolation 

Instead of direct statistical inference from datasets one can seek to reconstruct the 
underlying continuous field(s). For a meaningful analysis and interpretation of spa- 
tial data it is indeed often imperative and/or preferrable to use methods of parame- 
terization and interpolation to obtain estimates of the related field values throughout 
the sample volume. The reconstructed continuous field may subsequently be pro- 
cessed in order to yield a variety of interesting parameters. 

It involves issues of smoothing and spatial interpolation of the measured data 
over the sample volume, of considerable importance and interest in many different 
branches of science. Interpolation is fundamental to graphing, analysing and un- 
derstanding of sp atial da t a. Key referenc e s on t h e involved prob lems and solutions 
include those by iRipleyl (Il98ll) : IWatsonI d 1 992b : ICressid (Il993b . While of consid- 
erable importance for astronomical purposes, many available methods escaped at- 
tention. A s ystematic treatment and discussion within the astronomical context is 
the study bv lRvbicki & PressI (119921) . who focussed on linear systems as they devel- 
oped various statistical procedures related to linear prediction and optimal filtering, 
commonly known as Wiener filtering. An extensive, systematic and more general 
survey of available mathematical methods can be found in a set of publications by 
Lombardi & Schneider! (EoOlll2002lEo03l) . 
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Fig. 10. DTFE processed image of the Cosmic Web: GIF N-body simulation of structure 
formation in a ACDM cosmology. Part of the central X-slice. Simulation courtesy: J. Colberg. 



A particular class of spatial point distributions is the one in which the point pro- 
cess forms a representative reflection of an underlying smooth and continuous den- 
sity/intensity field. The spatial distribution of the points itself may then be used to 
infer the density field. This forms the basis for the interpretation and analysis of the 
large scale distribution of galaxies in galaxy redshift surveys. The number density 
of galaxies in redshift survey maps and N-body particles in computer simulations is 
supposed to be proportional to the underlying matter density. 



5.2 Local Interpolation: Natural Neighbour Methods 

The complex spatial geometry and large density variations marking the cosmic web 
ideally should be analyzed by a technique which would (1) not lose information 
against the backdrop of a highly inhomogeneous spatial resolution and (2) which 
is capable of tracing hierarchically structured and anisotropic spatial patterns in an 
entirely objective fashion. Nearly all existing techniques for analyzing galaxy red- 
shift surveys or numerical simulations of cosmic structure formation have important 
shortcomings with respect to how they treat the weblike geometry of the large scale 
matter distribution and trace the cosmic matter distribution over a wide range of 
densities. The limited available mathematical machinery has often been a major ob- 
stacle in exploiting the potentially large information content of the cosmic web. 

The various aspects characterizing the complex and nontrivial spatial structure 
of the cosmic web have proven to be notoriously difficult to quantify and describe. 
For the analysis of weblike patterns the toolbox of descriptive measures is still 
largely ad-hoc and is usually biased towards preconceived notions of their morphol- 
ogy and scale. None of the conventional, nor even specifically designed, measures 
of the spatial matter distribution have succeeded in describing all relevant features 
of the cosmic web. Even while they may succeed in quantifying one particular key 
aspect it usually excludes the ability to do so with other characteristics. 

For many applications a 'local' interpolation and reconstruction method appears 
to provide the prefeiTed path. In this case the value s of the variable at any point de- 



pend only on the data in its neighbourhood (see e.g lSambridge. Braun & McQueen , 



1995h . Local schemes usually involve a discretization of the region into an adaptive 
mesh. In data interpolation this usually represents a more realistic approach, and 
generically it also tends to be independent of specific model assumptions while they 
are very suited for numerical modeling applications. When the points have a regular 
distribution many local methods are available for smooth interpolation in multidi- 
mensional space. Smooth, local methods also exist for some specific irregular point 
distributions. A telling example are the 'locally constrained' point distributions em- 
ployed in applications of the finite element method. 

In this review we specifically concentrate on a wide class of tessellation-based 
multidimensional and entirely local interpolation procedures, commonly known 
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Fig. 11. Cosmic density field illustrating the large dynamic range which is present in the large 
scale matter distribution. In the left-hand frame the density field in a 10/r'Mpc wide slice 
through a cosmological simulation is depicted. In the subsequent frames zoom-ins focusing 
on a particular structure are shown. On all depicted scales structures are present. 



as Natural Neighbour Interpolatio n dWatsonl Il992t iBraun & Sambridgg, 1 1995 



Sukumai. 1998: Okabe 



mterpoiatio , 
et al.i l200d 



ch. 6). The local natural neighbour methods 
are based upon the Voronoi and Delaunay tessellations of the point sample, basic 



conce pts from the field of stochastic and computational geometry (see lOkabe et al 



200d and references therein). These spatial volume-covering divisions of space into 



mutually disjunct triangular (2-D) or tetrahedral (3-D) cells adapt to the local den- 
sity and the local geometry of the point distribution (see fig. [TO] and fig. fTTT i. The 
natural neighbour interpolation schemes exploit these virtues and thus adapt auto- 
matically and in an entirely natural fashion to changes in the density or the geom- 
etry of the distribution of sampling points. For the particular requirements posed 
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by astronomical and cosmological datasets, for which it is not uncommon to in- 
volve millions of points, we have developed a Unear first-order version of Nat- 
ural Nei ghbour Int erpolation, the Del aunay Tessellation Field Estimator ( DTFE, 
Bem ar deau & van d e Weygaert, 1996; Schaap & van de Wevga ert, 2000; Schaapl 
2007b . Instead of involving user-defined filters which are based upon artificial 
smoothing kernels the main virtue of natural neighbour methods is that they are 
intrinsically self-adaptive and involve filtering kernels which are defined by the lo- 
cal density and geometry of the point process or object distribution. 



5.3 Meshless methods 



The Natural Neighbour schemes, including DTFE, are mesh-based. With current 
technology this is computationally feasible as long as the domain of the matter 
density field has at most three dimensions. However, also interesting would be to 
extend attention to six-dimensional phase-space, incorporating not only the loca- 
tion of galaxies/particles but also their velocities. This will double the number of 
dimensions and makes mesh-based methods a real challenge. While the applica- 
tion of DTFE to the anal ysis of the p hase-space of da rk haloes has been shown to 
lead to very good results (lArad et al.L 12004) . studies by lAscasibar & Binneyl (l2005h 
and Sharma & Steinmetz ( 20061) argued it would be far more e fficient and rea- 
sonably accurate to resort to the simpler construct of a k-d tree (IBentlevl 1 1975 ; 
Friedman. Bentlev & Finkel , 1977 : Bentlev & Friedmanni 1978 : van de Wevgaert , 
19871 for an early astronomi cal implementation). While lXscasibar & Binnevl(l2005h : 
Sharma & Steinmetz (l2006h do not take into account that phase-space is not a sim- 
ple metric space but a symplectic one, it may indeed be a real challenge to develop 
mesh-based methods for the analysis of structure in phase-space. Although algo- 
rithms for the computation of higher dimensional Voronoi and Delaunay tessela- 
tions have been implemented (e.g., in CGAL ), the high running time and memory 
use make further processing computationally infeasible. 

Meshless spatial interpolation and approximation methods for datasets in spaces 
of dimensions greater than three may therefore provide the alternative of choice. 
There are a variety of meshl ess multivariate data interpolation schemes. Exam- 



ples are Shepard's interpolant (Shepardll 19681) . moving least squa res app roximants 



([Lancaster & Salkauskaslll981h . or Hardy's multiquadrics (IHardvtll971h . 



Spline interpolation 



Spline interpolation (Schoenberg 1946allbl) is based on interpolating between sam- 
pling points by means of higher-order polynomials. The coefficients of the poly- 
nomial are determined 'slightly' non-locally, such that a global smoothness in the 
interpolated function is guaranteed up to some order of derivative. The order of the 
interpolating polynomials is arbitrary, but in practice cubic splines are most widely 
used. Cubic splines produce an interpolated function that is continuous through the 
second derivative. To obtain a cubic spline interpolation for a dataset of 4- 1 points 
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Fig. 12. Interpolation of simulated Mars Orbiter Laster Altimeter (MOLA) data acquired 
from a digital elevation model of the Jokulsa and Fjollum region of Iceland, in anticipation of 
landscape reconstructions of terrain on planet Mars. The "original" highly resolved image of 
the terrain is shown in the top left frame. The comparison concerns data that were measured 
at the track points indicated in the top centre image. The medium resolution (209 x 492) 
interpolations are: natural neighbour (top right), linear (DTFE) interpolation (bottom left), 
nearest neighbour interpolation (bottom centre) and spline interpolation (bottom right). Evi- 
dently, the nn-neighbour map is the most natural looking reconstruction. Image courtesy of 
O. Abramov and A. McEwen, figure modified from Abramov & McEwen 2004 



separate cubics are needed. Each of these cubics should have each end point 
match up exactly with the end points to either side. At the location of these points 
the two adjacent cubics should also have e qual first and second deriva t ives. A full 
mathe matical derivation can be found in e.g. lGerald & Wheatlevi lll999h :l iPress et alj 
(Il992h . 

Spline interpolation is a widely used procedure. Equalising the derivatives has 
the effect of making the resulting interpolation appear smooth and visually pleasing. 
For this reason splines are for example frequently used in graphical applications. 
Splines can provide extremely accurate results when the original sample rate is no- 
tably greater than the frequency of fluctuation in the data. Splines however cannot 
deal very well with large gaps in the dataset. Because the gap between two points 
is represented by a cubic, these result in peaks or troughs in the interpolation. Also, 
splines are rather artificially defined constructs. 



Radial Basis Functions 



One of the most prom i sings schemes may be that of Rad ial Basis Functions, (RBFs, 
see e.g. IPowellL 1 1 9921: 1 Arnol'l 1 1 99 ll: l^ndlandl Eoosi) . RBFs may be used to de- 
termine a smooth density field interpolating three-dimensional spatial and four- 
dimensional spatio-temporal data sets, or even data sets in six-dimensional phase 
space. In the first step of this approach the implicit function is computed as a 
linear combination of translates of a single radial basis function. This function is 
determined by the geometric constraint that the input sample points belong to its 
zero set. If the input is a density map, the geometric constraint boils down to the 
implicit function interpolating the densities at the input points (and some addi- 
tional constraints preventing the construction of the zero function). The construc- 
tion of Radial Basis Functions with suitable interpolation properties is discussed in 
[Schab ack & Wendland (2001), while an early r eview o f the m athematical problems 
related to RBF-interpolation may be found in Powelll ( 1992 ). A nice overview of 
the state of the art in scattered data modeling using Radial Basis Functions may be 



obtained from the surveys iBuhmannI ( l2001l) : llske (2002); Lodha & Franke ( 1999). 

In practice variational implicit surfaces, based on Radial Basis Functions which 
mini mize some glo bal energy or cur vature functional, turn out to be very flex- 
ible dPinh. Turk & Slabaugh. .2002; .Turk & O'Brienl Il999h : they are adaptive to 
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curvature variations, can be used for enhancement of fine detail and sharp features 
that are missed or smoothed out by other implicit techniques, and can overcome 
noise in the input data since they are approximating rather than interpolating. Espe- 
cially the use of parameter dependent or anisotropic radial basis functions allows 
for g raceful treatment of sh arp features and provide multiple orders of smooth- 
ness dTurk&O'BrienLbOOlh . 



6 Spatial Tessellations 

6.1 Stochastic and Computational Geometry 

Random spatial tessellations are a fundamental concept in the fields of Stochastic 
Geometry and Computational Geometry. 

Stochastic Geometry, or geometric probabiUty theory, is the subject in math- 
ematics concerned with the problems that arise when we ascribe probability dis- 
tributions to geometric objects such as points, lines, and planes (usually in Eu- 
clid ian spaces), or to geometric operations such as rotations or projections (see 
e.g. Stovan. Kendall & Mecke , 19871 1995b. A formal and restricte d definition of 



stochastic geometry was given by IStovan. Kendall & Meckel (119870 . who defined 



the field as the branch of mathematics devoted to the study of geometrical struc- 
tures whcih can be described by random sets, in particular by point processes, in 
suitable spaces. Since many problems in stochastic geometry are either not solvable 
analytically, or only in a few non-generic cases, we have to resort to the help of the 
computer to find a (partial) answer to the problems under consideration. This makes 
it necessary to find efiicient algorithms for constructing the geometrical objects in- 
volved. 

Computational Geometry is the branch of computer science that is concerned 
with finding the computational procedures for solving geometric problems, not just 

the geo metric problems arising from stochastic geometry but geometric problems in 

gener al jBoissonnat. Yvinec. Bronnmimanill998l:lde Berg et al. [ |2000l: l iGoodman & O'Rourkel 



20041) . It is concerned with the design and analysis of algorithms and software for 
processing geometric objects and data. Typical problems are the construction of spa- 
tial tessellations, like Voronoi diagrams and Delaunay meshes, the reconstruction of 
objects from finite point samples, or finding nearest neighbours in point sets. Meth- 
ods from this field have many applications in applied areas like Computer Graphics, 
Computer Vision, Robotics, Computer Aided Design and Manifacturing. 



6.2 Random Tessellations 

Random tessellations or mosaics occur as primary objects of various processes of 
division of space ^R'' into convex cells or as secondary or auxiliary objects in e.g. 
various statistical problems. 

Simply stated, a tessellation is an arrangement of polytopes (2-D: polygons, 
3-D: polyehdra) fitting together without overlapping so as to cover d-dimensional 
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Fig. 13. A full 3-D tessellation comprising 1000 Voronoi cells/polyhedra generated by 1000 
Poissonian distributed nuclei. Courtesy: Jacco Bankers 
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Fig. 14. An illustration of a distribution of nuclei (stars) in a square (top) and its correspond- 
ing Delaunay triangulation (bottom left) and Voronoi tessellation (bottom right), assuming 
periodic boundary conditions. 

space ^R'' {d - 1 , 2, . . .), or a subset X c !R''. Usually one requires that the cells are 
convex and compact with disjoint interiors, although tessellations may also involve 
non-convex cells (an example are Johnson-Mehl tessellations). Posing the additional 
requirement of convexity implies that all interfaces separating pairs of cells are pla- 
nar (for 3-D), and all edges common to the boundaries of three or more cells are 
linear, implying each cell to be a (convex) polyhedron. Formally, a tessellation in 
^R'' is a set 7" = {X,) of c/-dimensional sets X, c ^R'' called cells such that 
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Xi n X; = 



U 



J - " for / 5t 



(1) 



#{Xi € r ; X; n B ?i 0) < oo V bounded B c 



with iti the interior of cell X, ( lM0lleii Il989i Il994 . The first property implies the 
interiors of the cells to be disjoint, the second one that the cell aggregate {Xi} is 
space-filling, and the third one that 7~ is a countable set of cells. 



6.3 Voronoi Tessellation 

The Voronoi tessellation of a point set !P is the division of space into mutually 
disjunct polyhedra, each Voronoi polyhedron consisting of the part of space closer to 
the defining point than any of the other points (iVoronoL.1908; Oka be et al., 2000i). 
Assume that we have a distribution of a countable set P of nuclei {x,) in ^R''. Let 
Xi , X2, X3, . . . be the coordinates of the nuclei. Then the Voronoi region 'V, of nucleus 
i is defined by the points x (fig.EJi, 



"Vi = {x\d{x, Xi) < d{x, xj) Vj ^ /} 



(2) 



where d(x, y) is the Euclidian distance between x and y. In other words, is the set 
of points which is nearer to x,- than to Xj, j + i. From this basic definition, we can 
directly infer that each Voronoi region 'V, is the intersection of the open half-spaces 
bounded by the perpendicular bisectors (bisecting planes in 3-D) of the line seg- 
ments joining the nucleus / and any of the the other nuclei. This implies a Voronoi 
region 'Vi to be a convex polyhedron (a polygon when in 2-D), a Voronoi polyhe- 
dron. Evidently, the concept can be extended to any arbitrary distance measure (Icke, 
priv. comm.). The relation between the point distribution f and its Voronoi tessella- 
tion can be clearly appreciated from the two-dimensional illustration in fig. [14] The 
complete set of Voronoi polyhedra constitute a space-filling tessellation of mutually 
disjunct cells, the Voronoi tessellation 'V{P) relative to !P. A good impression of the 
morphology of a complete Voronoi tessellation can be seen in fig. [13] a tessellation 
of 1000 cells generated by a Poisson distribution of 1000 nuclei in a cubic box. The 
Voronoi foam forms a packing of Voronoi cells, each cell being a convex polyhedron 
enclosed by the bisecting planes between the nuclei and their natural neighbours. 
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3-D Voronoi Tessellation Elements 

(see also fig.ll5t 





Voronoi cell 


- Polyhedron 

- defined by nucleus / e f 

- volume of space closer to i 
than any other nucleus meV 




Voronoi wall (facet) 


- Polygon 

- defined by nuclei j) e P 

- all points x with equal distance to j) 

& larger distance to any other nucleus m e f 

- constitutes part surface cells: "Vj, "Vj 


Aijk 


Voronoi edge 


- Line segment 

- defined by nuclei j, k) eV 

- all points x with equal distance to (;, j, k) 

& larger distance to any other nucleus m e f 

- constitutes part rim Voronoi cells: 'Vj, "Vj, "Vk 
constitutes part rim Voronoi walls: Zjj, and Zj^ 




Voronoi vertex 


- Point 

- defined by nuclei j, k,l) eP 

- equidistant to nuclei (i, j, k, I) 

- closer to (/, j, k, I) than to any other nucleus m eV 

- circumcentre of (Delaunay) tetrahedron /, k, I) 



Voronoi elements 

Taking the three-dimensional tessellation as the archetypical representation of struc- 
tures in the physical world, the Voronoi tessellation 'V{'P consists of four constituent 
elements: Voronoi cells, Voronoi walls, Voronoi edges and Voronoi vertices. Ta- 
ble |63] provides a listing of these elements together with a description of their re- 
lation with respect to the nuclei of the generating point set V, augmented by the 
accompanying illustration in fig. [15] 



Generalized Voronoi Tessellations 



The Voronoi tessellation can be generalized. iMiles d 19701) defined the generalized 
Voronoi tessellation 'V„. The original Voronoi tessellation is 'V\ — 'V. 

This involves the extension of the definition of the Voronoi cell 'V, generated by 
one nucleus / to that of a higher-order Voronoi cell "V^X/i, . . . , ik) generated by a set 
of k nuclei {i\,. . . , ik] e P. Each k - order Voronoi cell 'V^{i\, . . . , ik) consists of 
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Fig. 15. The four Voronoi elements of a Voronoi tessellation generated by a nucleus set P. 
See tablelO 



that part of space in which the p oints x have the k nuclei { i\ , . . . , ik) e !P as their k 
nearest neighbours. In addition to iMilesI (fl970[ 1 1 972l 1 1 982h see lOkabe etal] (l2000i 
chapter 3) for more details and references. 



Voronoi Uniqueness 

An inverse look at tessellations reveals the special and degenerate nature of Voronoi 
tessellations. Given a particular tessellation one might wonder whether there is a 
point process which would have the tessellation as its Voronoi tessellation. One 
may demonstrate that in general this is not true. By defining a geometric procedure 
to reconstruct the generating nucleus distribution one may straightforwardly infer 
that there is not a unique solution for an arbitrary tessellation. 
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This may be inferred from the study by Chiu. van de Wevgaert & Stovan ( 19961) . 
who defined and developed a nucleus reconstruction procedure to demonstrate that a 
two-dimensional section through a three- or higher-dimensional Voronoi tessellation 
will itself not be a Voronoi tessellation. By doing so their work clarified the special 
and degenerate nature of these tessellations. 



6.4 the Delaunay Tessellation 

Pursuing our census of Voronoi tessellation elements ftable l6.3l l. we found that each 
set of nuclei /, j, k, I corresponding to a Voronoi vert ex £)(/, i,k , l) defi nes a unique 
tetrahedron. This is known as Delaunay tetrahedron ( lDelaunavill934 ). 



Each Delaunay tetrahedron is defined by the set of four points whose circum- 
scribinn sphere d oes not contain any of the other points in the generating set 
dDelaunavi Il934i triangles in 2D: see fig. [T4l i. For the countable set P of points 
{x,} in 51'', a Delaunay tetrahedron £),„ is the simplex T defined by (1 -i- d) points 
{x,i, . . . ,x,y+i)) e V (the vertices of this cZ-dimensional tetrahedron) such that the 
corresponding circumscribing sphere Smiym) with circumcenter C,„ and radius 
does not contain any other point of V, 



'Dm - T(Xii, . . . , X/(^+i)) 


with d(Cm,Xj) 


> Rm 




V j^il,... 


i{d + 1) 



(3) 



Following this definition, the Delaunay tessellation of a point set f is the uniquely 
defined and volume-covering tessellation of mutually disjunct Delaunay tetrahedra. 

Figure [14] depicts the Delaunay tessellation resulting from a given 2-D distribu- 
tion of nuclei. On the basis of the figure we can immediately observe the intimate 
relation between a Voronoi tessellation and its Delaunay tessellation. The Delaunay 
and Voronoi tessellations are hke the opposite sides of the same coin, they are each 
others dual: one may directly infer one from the other and vice versa. The combina- 
torial structure of either tessellation is completely determined from its dual. 

The duality between Delaunay and Voronoi tessellations may be best appreci- 
ated on behalf of the following properties: 

• Circumcentre Gr Voronoi vertex: 

The center of the circumsphere of a Delaunay tetrahedron is a vertex of the 
Voronoi tessellation. This follows from the definition of the Voronoi tessellation, 
wherein the four nuclei which form the Delaunay tetrahedron are equidistant 
from the vertex. 
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• Contiguity condition: 

The circumsphere of a Delaunay tetrahedron is empty and cannot contain any 
nucleus in the set !P. If there would be such a nucleus it would be closer to the 
centre than the four tetrahedron defining nuclei. This would render it impossible 
for it being the vertex of all corresponding Voronoi cells. 



Natural Neighbours 

A pair of nuclei / and j whose Voronoi polyhedra and 'Vj have a face in common 
is called a contiguous pair and a member of the pair is said to be contiguous to the 
other member. Contiguous pairs of nuclei are each other's natural neighbour. 

Natural neighbours of a point / are the points j with whose Voronoi 'Vj its 
Voronoi cell 'V, shares a face or, in other words the points with which it is con- 
nected via a Delaunay tetrahedron. This unique set of neighbouring points defines 
the neighbourhood of the point and represents the cleanest definition of the sur- 
roundings of a point (see fig.[T6]l, an aspect which turns out to be of seminal impor- 
tance for the local interpolation method(s) discussed in this contribution.. 




Fig. 16. The dual relationship between Voronoi (solid) and Delaunay (dashed) tessellations of 
a set of nuclei (circles). Left: Zoom-in on the three Delaunay triangles corresponding to a set 
of five nuceli (black dots) and the corresponding Voronoi edges. The circle is the circumcircle 
of the lower Delaunay triangle, its centre (an^ow) is a vertex of the Voronoi cell. Note that the 
circle does not contain any other nucleus in its interior ! Right: a zoom in on the Voronoi cell 
•y, of nucleus / (black dot). The Voronoi cell is surrounded by its related Delaunay triangles, 
and clearly delineates its Natural Neighbours (open circles). 
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6.5 Voronoi and Delaunay Statistics 



In particular for practical applications the knowledge of statistical properties of 
Voronoi and Delaunay tessellations as a function of the generating stochastic point 
processes is of considerable interest. However, despite their seemingly simple def- 
inition it has proven remarkably difficult to obtain solid analytical expressions for 
statistical properties of Voronoi tessellations. Moreover, with the exception of some 
general tessellation properties nearly all analytical work on Voronoi and Delaunay 
tessellations has concentrated on those generated by homogeneous Poissonian point 
distributions. 

Statistical knowledge of Delaunay tessellations generated by Poissonian nuclei 
in d- dimensional space is relatively complet e. Som e important d i stribu ti on func- 
tions are known, mainly du e to the work of Miles ( 1970i 1972 , 1974); Kendall 
19991) : M0llei ( 1989 , 1994 ). For Voronoi tessellations, even for Poissonian nuclei 



analytical results are quite rare. Only very few distribution functions are known, 
most results are Umited to a few statistical moments: expectation values, vari- 
ances an d correlation coefficients. Most of t hese results s te m frorn the pi o neering 
works of lMeijerind (Il953k : lOilbertl (Il962l) : iMilesI (Il970l) : lM0lleil (Il989l) . lM0ller 
(Il989h provides analytical formulae for a large number of first-order moments 
of c/-dimensional Poisson- Voronoi tessellations, as well as of i-dimensional sec 



tion through them, fol lowing the work of iMilesI (119721 11974 , IT982D (also see 
van de Weygaertl 1991). Numerical Monte Carlo evaluations have proven te be the 
main source for a large range of addi tional statistical resu lts. For an extensive listing 
and discussion we refer the reader to Okabe et al. I (l2000l) . 

Of prime importance for the tessellation interpolation techniques discussed in 
this review are two fundamental characteristics: (1) the number of natural neigh- 
bours of the generating nuclei and (2) the volume of Voronoi and Delaunay cells. 

In two dimensions each point h as on average exactly 6 natural neighbours (see 
e.g. llcke & van de Wevgaertlll987b . irrespective of the character of the spatial point 
distribution. This also implies that each point belongs to on average 6 Delaunay tri- 
angles. Going from two to three dimensions, the character of the tessellation changes 
fundamentally: 

• Dependence Point Process: 

the average number of natural neighbours is no longer independent of the un- 
derlying point distribution 

• Integer. 

The average number of natural neighbours is not an integer: 
for a Poisson distribution it is ~ 13.4 ! 

• Delaunay tetrahedra/Voronoi Vertices: 

For a Poisson distribution, the # vertices per Voronoi cell ~ 27.07. 
While in the two-dimensional case the number of vertices per Voronoi cell has 
same value as the number of natural neighbours, in three dimensions it is en- 
tirely different ! 
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As yet it has not been possible to derive, from first principle, a closed analytical 
expression for the distribution functions of the volumes of Voronoi polyhedra and 
Delaunay cells in a Poisson Voronoi tessellation. However, the fitting formula sug- 
gested by Kiang ( 1966) has proven to represent a reasonably good approximation. 
Accordingly, the probability distribution of the volume of of a Voronoi polyhedron 
in -dimensional space follows a gamma distribution 



1: 



r{q) Viv-v) 



exp -q 



(4) 



w ith Vy the s ize of the Voronoi cell and (V^) the average cell size. The conjecture of 
iand (119661) is that the index has a value q - 2d for a tessellation in tZ-dimensional 
space (ie. q = A for 2-D space and q - b for 3-D space). Even though other studies 
indicated slightly different valu es, we have fou nd the Kiang suggestion to be quite 
accurate (see sect. ll 1.2l also see lSchaapl (l2007h l 

While the distribution for Voronoi cell volumes involves a conjecture, for De- 
launay cells D it is possible to derive the ergodic distribution from first principle. Its 
d + \ vertices completely specify the Delaunay tetrahedron £). If c, R its circumra- 
dius, its vertices are the points {c-i-/?u,). In this the unit vectors (u,) (/ = 0, . . . , cZ-i- 1), 
directed towards the vertices, determine the shape of the Delaunay tetrahedron. 
Given that is the volume of the unity simplex {hq, . . . , hj), and the volume Vj) of 
the Delaunay tetrahedron given by 



(5) 



MilesI d 1 974 : lM0lleil ( 1 1 989l 1994b found that for a Poisson point process of intensity 
n in ^/-dimensional space the distribution is specified by 



fM - M{uo,...,Ud),R) 



= a(n,d)Ad R'^''^ exp(-«w,;/?'') , 



(6) 



with a>d the volume of the unity sphere in c/-dimensional space (ai2 - 2k and 
cji = An) and a{n, d) is a constant dependent on number density n and dimension d. 
In other words, the circumradius R of the Delaunay tetrahedron is ergodically inde- 
pendent of its shape, encapsulated in A d. From this one finds that the the distribution 
law of iicodR^' is thex2d/'^ distribution ([Kendall, Il999h : 



f(R)dR 



(noJdR'y 



'd\d-l 



exp {-najdR"^) dinajdR"^) . 



(7) 



{d-iy. 

It is of considerable importance to realize that even for a uniform density field, 

represented by a discrete point distribution of intensity n, neither the Voronoi nor the 
Delaunay tessellation will involve a volume distribution marked by a considerable 
spread and skewness. 
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Fig. 17. The Delaunay tessellation of a point distribution in and around a filamentary feature. 
The generated tessellations is shown at three successive zoom-ins. The frames form a testi- 
mony of the strong adaptivity of the Delaunay tessellations to local density and geometry of 
the spatial point distribution. 

6.6 Spatial Adaptivity 

The distribution of the Delaunay and Voronoi cells adjusts itself to the characteris- 
tics of the point distribution: in sparsely sampled regions the distance between the 
natural neighbours is large, and also if the spatial distribution is anisotropic this 
will be reflected in their distribution. This may be readily appreciated from fig.fTTl 
At three successive spatial scales the Delaunay tessellations has traced the density 
and geometry of the local point distribution to a remarkable degree. On the basis of 
these observations it is straightforward to appreciate that the corresponding Delau- 
nay tessellation forms an ideal adaptive multi-dimensional interpolation grid. 
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Note that not only the size, but also the shape of the Delaunay simplices is fully 
determined by the spatial point distribution. The density and geometry of the lo- 
cal point distribution will therefore dictate the resolution of spatial interpolation and 
reconstruction procedures exploiting the Delaunay tessellation. The prime represen- 
tatives of such methods are the Natural Neighbour techniques (see sect. |7]) and the 
Delaunay Tessellation Field Estimator. 

These techniques exploit the fact that a Delaunay tetrahedron may be regarded 
as the optimal multidimensional interpolation interval. The corresponding minimal 
coverage characteristics of the Delaunay tessellation thus imply it to be optimal for 
defining a network of multidimensional interpolation intervals. The resulting inter- 
polation kernels of Natural Neighbour and DTFE interpolation not only embody an 
optimal spatial resolution, but also involve a high level of adaptivity to the local 
geometry of the point distribution (see sect. |7] and sect. [8]l. Within this context it is 
no coincidence that in many computer visualisation applications Delaunay tessella- 
tions have acquired the status of optimal triangulation. Moreover, the superb spatial 
adaptivity of the volumes of Delaunay and Voronoi polyhedra to the local point den- 
sity may be readily translated into measures for the value of the local density. This 
forms a crucial ingredient of the DTFE formalism (see sect. [8]). 



6.7 Voronoi and Delaunay Tessellations: context 



T he earliest signifi cant use of Vo r onoi r egions seems to have occurred in the work 
of iDirichleli dlSSOl) and IVoronoil (Il908h in their investigations on the reducibility 
of positive definite quadratic forms. However, Dirichlet and Voronoi tessellations 
as applied to random point patter ns appear to have ar isen independently in various 



fields of science and technology dOkabe et al. . I2OOOI) . For example, in crystallog 



raphy, one simple model of crystal growth starts with a fixed collection of sites in 
two- and three-dimensional space, and allows crystals to begin growing from each 
site, spreading out at a uniform rate in all directions, until all space is filled. The 
"crystals" then consist of all points nearest to a particular site, and consequently are 
just the Voronoi regions for the original set of points. Similarly, the statistical analy- 
sis of metereologic al data led to the formulation of Voronoi regions under the name 
Thiessen polygons (lThiessenlll91 ll) . 

Applications of Vorono i tessellations can there fore be found in fie lds as diverse 

as agriculture and forest ry (Fisch er & Milesl 1973b. astrophysics (e.g. Kiand, 1966 ; 

Icke & van de Wevgaerttil987MEbeUng & Wiedenmanni 19931: Bernardeau & van de Wevgaert , 



1996t Molchanov. Surgailis & Wovczvnskil 1997 : Schaap & van de Wevgaertll200C : 
Cappellari & CopinL 2003 Ritzerveld & IckeL l2006[). ecologv. zoo lo gy and botany. 



cell bio logy ( Gor et al. , 20061). protein research jLiang et al.L 1998allbl: Liang 



1998d) . cancer research (IKansal et al. , I2OOOI: ISchaller & Mever-HermannI 12005^ 
chemistry, crystal growth and struct ure, materials science (see e.g. Torguato 12002 , 
incl. many references), geophysics dSambridge. Braun & McQueen , 1995), geog- 
raphy and geographic informat ion systems ( Boots , 19841 : Gold. Remmele & Roos , 
19971). communication theory ( Baccelli & Zuyev 19991) and art and archaeology 
dKimia & Leymariell2001uLey marie & Kimia I2003ir ' )ue to the diversity of these 



Ldelsbrunner & Woodward . 
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applications it has acquired a set of alternative names, such as Dirichlet regions, 
Wigner-Seitz cells, and Thiessen figures. 



7 Natural Neighbour Interpolation 



The Natural Neighbour Interpolation formalism is a generic higher-order multidi- 
mensional interpolation, smoothing and modelling procedure utilizing the concept 
of natural neighbours to obtain locally optimized measures of system char acteris- 
tics. Its theoretical basis was developed and introduced by ISibson (Il98l[). while 
extens i ve treatments and elaborations of nn-interpolation may be found in IWatson 
(ll992h : ISukumarl ( Il998h . Natural neighbour interpolation produces a conservative. 



2 
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Fig. 18. Natural Neighbour Interpolants. Example of NN interpolation in two dimensions. 
Left: the Voronoi cell 'Vi's. generated by a point x. Right: the 2nd order Voronoi cell 'V2('s., Xi ), 
the region of space for which the points x and xj are the closest points. Image courtesy M. 
Sambridge and N. Sukumar. Also see Sambridge, Braun & McQueen 1995 and Sukumar 
1998. 



artifice-free, result by finding weighted averages, at each interpolation point, of the 
functional values associated with that subset of data which are natural nei ghbors of 



each interpolation point. Unlike other schemes, like Shepard's interpolant jShepard , 



19681) . where distance-based weights are used, the Sibson natural neighbour inter- 
polation uses area-based weights. According to the nn-interpolation scheme the in- 
terpolated value /(x) at a position x is given by 

/(X) = 2 0„„,,(x)/;-, (8) 
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in which the summation is over the natural neighbours of the point x, i.e. the sample 




Fig. 19. Natural Neighbour Interpolants. Example of NN interpolation in two dimensions. 
The basis of the nn-interpolation kernel i^,,,, i , the contribution of nucleus Xj to the interpolated 
field values. Image courtesy M. Sambridge and N. Sukumar. Also see Sambridge, Braun & 
McQueen 1995 and Sukumar 1998. 

points j with whom the order-2 Voronoi cells ^2(x,Xj) are not empty (fig. [T8l[T9] l. 
Sibson interpolation is based upon the interpolation kernel (p{x, Xj) to be equal to the 




Fig. 20. Natural Neighbour Interpolants. Example of NN interpolation in two dimensions. 
3-D illustration of the nn-interpolation kernel 0„„_i, the contribution of nucleus Xi to the in- 
terpolated field values. Image courtesy M. Sambridge, also see Braun & Sambridge 1995. 
Reproduced by permission of Nature. 
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Fig. 21. Natural Neighbour Kernels: illustration of the locally defined natural neighbour ker- 
nels 0„„j for four different nuclei j. Image courtesy M. Sambridge, also see Braun & Sam- 
bridge 1995. Reproduced by permission of Nature. 



normalized order-2 Voronoi cell, 



,,,(x) = 



J?12(X,X,) 



(9) 



in which ^(x) = 2y.^(x,xy) is the volume of the potential Voronoi cell of point x 
if it had been added to the point sample P and the volume J?l2(x, x,) concerns the 
order-2 Voronoi cell 'ViCx, x,), the region of space for which the points x and Xj are 
the closest points. Notice that the interpolation kernels <p are always positive and 
sum to one. 

The resulting function is continuous everywhere within the convex hull of the 
data, and has a continuous slope everywhere except at the data themselves (fig.l20l 
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fig.lSTTi. Beautiful two-dimensional examples of nn-interpolation applications testify 
of its virtues (see fig.[T2li. 

An interesting study of its performance, in comparison with other interpola- 
tion/approximation methods, concerned the test involving data acquired by the Mars 
Orbiter Laser Altimeter (MOLA), one of the instruments on the Mars Global Sur- 
veyor (deceased November 2006). Applying various schemes towards the recon- 
struction of a terrain near the Korolev crater, a large crater in the north polar region 
of Mars, demonstrated that nn-neighbou r interpolation produced the most impres- 
sive reproduction of the original terrain dAbramov & McEwenll2004l) . In compari- 
son to the other methods - including nearest neighbour intp., spline intp., and linear 
DTFE type interpolation - the nn-neighbour map not only looks most natural but 
also proves to contain fewer artefacts, both in number and severity. The test recon- 
structions of the Jokulsa and Fjollum region of Iceland in fig. [12] provide ample 
proof. 

Natural neighbour interpolation may rightfully be regarded as the most general 
and robust method for multidimensional interpolation available to date. This smooth 
and local spatial interpolation technique has indeed gradually acquired recognition 
as a dependable, optimal, local method. For the two-dimensional applications it 
has seen highly interesting appUcations in geophysical fluid dynam ics calculations 
dBraun & Sambridgei 1 19951 : ISambridge. Braun & McOueenl 1 1995b . and in equiva- 
lent schemes for solid mechanics problems ( Sukumai , 19981) . Both applications used 
n-n interpolation to solve partial differential equations, showing its great potential 
in the field of computational fluid mechanics. 

While almost optimal in the quality of its reconstructions, it still involves a heavy 
computational effort. This is certainly true for three or higher dimensional spaces. 
This has prodded us to define a related local, adaptive method which can be applied 
to large astronomical datasets. 



8 DTFE: the Delaunay Tessellation Field Estimator 

For three-dimensional samples with large number of points, akin to those found in 
large cosmological computer simulations, the more complex geometric operations 
involved in the pure nn-neighbour interpolation still represent a computationally 
challenging task. To deal with the large point samples consisting of hundreds of 
thousands to several millions of points we chose to follow a related nn-neigbhour 
based technique that restricts itself to pure linear int erpolation. DTFE u ses to ob- 
tain optimal local estimates of the spatial density (see Okabe et al. . 200d sect 8.5), 



while the tetrahedra of its dual Delaunay tessellation are used as multidimensional 
intervals for linear interp olation of the field v alues sampled or estimated at the loca- 
tion of the sample points ( lOkabe et al. [ |2000[ ch. 6). The DTFE technique allows us 



to follow the same geometrical and structural adaptive properties of the higher or- 
der nn-neighbour methods while allowing the analysis of truely large data sets. The 
presented tests in this review will demonstrate that DTFE indeed is able to highlight 
and analyze essential elements of the cosmic matter distribution. 
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8.1 DTFE: the Point Sample 

The DTFE procedure, outlined in the flow diagram in fig. |22] involves a sequence 
of steps. It starts with outlining the discrete point sample T in tZ-dimensional space 

{xi,...,x^), (10) 

In the applications described in this study we restrict ourselves to Euclidian spaces, 
in particular 2- or 3-dimensional Euclidian spaces. At the locations of the points 
in the countable set P the field values {/(x,), / = 1, . . . , A^) are sampled, or can be 
estimated on the basis of the spatial point distribution itself. The prime example of 
the latter is when the field / involves the density field itself. 

On the basis of the character of the sampled fi eld /(x) we need to distinguish 



two op tions. The first option is the one defined in iBernardeau & van de Wevgaerl 



(119961) . the straightforward multidimensional li near interpolation of measured field 



val ues on th e basis of the Delaunay tessellation. ISchaap & van de Wevgaertl ( l2000l) 



and Schaapl ( 2007.) extended this to the recovery of the density or intensity field. 



DTFE is therefore characterized by a second option, the ability to reconstruct the 

underlying density field from the discrete point process itself. 

The essential extension of DTFE wrt. iBemardeau & van de Wevgaerl (Il99d) IS 



that it allows the option of using the point sample process V itself as a measure for 
the value of the density at its position. The latter poses some stringent requirements 
on the sampling process. It is crucial for the discrete spatial point distribution to 
constitute a fair sample of the underlying continuous density field. In other words, 
the discrete point sample f needs to constitute a general Poisson point process of 
the density field. 

Such stringent requirements on the spatial point process V are not necessary 
when the sampled field has a more generic character As long as the sample 
points are spread throughout most of the sample volume the interpolation pro- 
cedure will yield a properly representative field reconstruction. It was this sit- 
uation with respect to cosmic velocity fields which lead to the original defini- 
tion of the Delaunay spatial interpolation procedure forming the basis of DTFE 
(IBernardeau & van de Wevgaertlll996t l IBernardeau et al.Lll997h . 



8.2 DTFE: Linear Interpolation 

At the core of the DTFE procedure is the use of the Delaunay tessellation of the dis- 
crete point sample process (see sect. 16.4b as an adaptive spatial linear interpolation 
grid. Once the Delaunay tessellation of V is determined, it is used as a multidimen- 
sional linear interpolation grid for a field /(r). That is, each Delaunay tetrahedron 
is supposed to be a region with a constant field gradient, V/. 

The linear interpolation scheme of DTFE exploits the same spatially adaptive 
characteristics of the Delaunay tessellation generated by the point sample P as that 
of regular natural neighbour schemes (see sect.|7] eqn.|8]i. For DTFE the interpola- 
tion kernel (pdt.A^) is that of regular linear interpolation within the Delaunay tetra- 
hedron in which x is located, 
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(11) 



in which the sum is over the four sample points defining the Delaunay tetrahedron. 
Note that not only the size, but also the shape of the Delaunay simplices is fully 
determined by the spatial point distribution. As a result the resolution of the DTFE 
procedure depends on both the density and geometry of the local point distribution. 
Not only does the DTFE kernel embody an optimal spatial resolution, it also in- 
volves a high level of adaptivity to the local geometry of the point distribution (see 
sect. lU.ll i. 

Also note that for both the nn-interpolation as well as for the linear DTFE inter- 
polation, the interpolation kernels 0, are unity at sample point location x, and equal 
to zero at the location of the other sample points j (e.g. see fig.|20l). 



if j. 



(12) 



where Xj is the location of sample point j. 

In practice, it is convenient to replace eqn.[TT]with its equivalent expression in 
terms of the (linear) gradient V/|^^ inside the Delaunay simplex m. 



/(X) = /(X,.) + V/ -(x-x,). 



(13) 



The value of V/|^^ can be easily and uniquely determined from the (I + D) field val- 
ues fj at the sample points constituting the vertices of a Delaunay simplex. Given 
the location r - {x,y,z) of the four points forming the Delaunay tetrahedra's ver- 
tices, ro, ri, r2 and rs, and the value of the sampled field at each of these locations, 
/o, /i, /a and fi, and defining the quantities 



Ax/j — xq , 

Ay„ ^ jn- yo ; 

/iz„ = Zn - ZO 



for « = 1,2,3 



(14) 



as well as Af,, = fn-fo (« - 1, 2, 3) the gradient V/ follows from the inversion 



V/ 



(15) 



dx 

df 
dy 

^dz' 

Once the value of Vy has been determined for each Delaunay tetrahedron in the 
tessellation, it is straightforward to determine the DTFE field value /(x) for any lo- 
cation X by means of straightforward linear interpolation within the Delaunay tetra- 
hedron in which x is located (eqn.flj]). 









'Axi Ay I Az\ 




Af2 


; A = 


Ax2 Ay 2 Az2 








,Ax3 Ay 3 Az-i, 
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The one remaining complication is to locate the Delaunay tetrahedron £)„, in 
which a particular point x is located. This is not as trivial as one might naively 
think. It not necessarily concerns a tetrahedron of which the nearest nucleu s is a ver- 
tex. Fortunately, a v ery efficient method, the walking triangle algorithm ( ILawson , 
1977t ISloanLll987h has been developed. Details of t he method may be found in 
Sambridge. Braun & McOueenI (ll995i : ISchaapl (l2007h . 



8.3 DTFE: Extra Features 



While DTFE in essence is a first order version of Natural Neighbour Interpolation 
procedure, following the same adaptive multidimensional interpolation characteris- 
tics of the Delaunay grid as the higher-order nn-neighbour techniques, it also incor- 
porates significant extensions and additional aspects. In particular, DTFE involves 
two extra and unique features which are of crucial importance for the intended cos- 
mological context and application: 

• Volume-weighted: 

The interpolated quantities are volume-weighted, instead of the implicit mass- 
weighted averages yielded by conventional grid interpolations. 

• Density estimates: 

The spatial adaptivity of the Delaunay/Voronoi tessellation to the underlying 
point distribution is used to estimate the local density. 



9 DTFE: the field reconstruction procedure 

The complete DTFE reconstruction procedure, essential steps of which are illus- 
trated in Fig. |23] can be summarized in terms of the flow diagram in Fig. l22]and 
consists of the foUowing sequence of steps. 

• Point sample 

Defining the spatial distribution of the point sample: 
+ Density field: 

point sample needs to be a general Poisson process of the (supposed) under- 
lying density field, i.e. it needs to be an unbiased sample of the underlying 
density field. 
+ General (non-density) field: 

no stringent requirements upon the stochastic representativeness of the sam- 
pling process will be necessary except that the sample volume is adequately 
covered. 

• Boundary Conditions 

An important issue, wrt the subsequent Delaunay tessellation computation and 
the self-consistency of the DTFE density and velocity field reconstructions, is 
that of the assumed boundary conditions. These will determine the Delaunay 
and Voronoi cells that overlap the boundary of the sample volume. Dependent 
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Fig. 22. Flow diagram illustrating the essential ingredients of the DTFE procedure. 
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Fig. 23. Summary: overview of the essential steps of the DTFE reconstruction procedure. 
Given a point distribution (top left), one has to construct its corresponding Delaunay tessel- 
lation (top right), estimate the density at the position of the sampling points by taking the 
inverse of the area of their corresponding contiguous Voronoi cells (bottom right) and fi- 
nally to assume that the density varies linearly within each Delaunay triangle, resulting in a 
volume-covering continuous density field (bottom left). 




upon the sample at hand, a variety of options exists: 

+ Vacuum boundary conditions: 

outside the sample volume there are no points. This will lead to infinitely 
extended (contiguous) Voronoi cells surrouding sample points near the 
boundary. Evidently, these cells cannot be used for DTFE density field es- 
timates and field interpolations: the volume of the DTFE reconstruction is 



52 



Rien van de Weygaert & Willem Schaap 



smaller than that of the sample volume. 



+ Periodic boundary conditions: 

the point sample is supposed to be repeated periodically in boundary boxes, 
defining a toroidal topology for the sample volume. The resulting Delaunay 
and Voronoi tessellations are also periodic, their total volume exactly equal 
to the sampel volume. While specific periodic tessellation algorithms do 
exist (Ivan de Weygaern . 1 19941) . this is not yet true for most available rou- 
tines in standard libraries. For the analysis of N-body simulations this is the 
most straightforward and preferrable choice. 



+ Buffer conditions: 

the sample volume box is surrounded by a bufi'erzone filled with a syn- 
thetic point sample. The density of the synthetic buffer point sample should 
be related to the density in the nearby sample volume. The depth of the 
bufferzone depends on the density of the synthetic point sample, it should 
be sufficiently wide for any Delaunay or Voronoi cell related to a sample 
point not to exceed the bufferzon e. A clear condition for a sufficiently d eep 
bufferzone has been specified by Neyrinck. Gnedin & Hamiltonl ( 2005 ). 



When involving velocity field analysis, the velocities of the buffer points 
should also follow the velocity distribution of the sample and be in accor- 
dance with the continuity equation. Relevant examples of possible choices 
are: 



- internal: the analyzed sample is a subsample embedded within a large 
sample volume, a sufficient number of these sample points outside the anal- 
ysis volume is taken along in the DTFE reconstruction. 



- rand om cloning technique: akin to the technique described by lYahil et al 
(Il99lh . 



- constrained r andom field: realizations employing the existing correla- 

tions in the field (lBertschingeiill987l:lHoff^man & Ribaklll99ll: Ivan de Weygaert & Bertschingen, 



1996HZaroubi et al.Lll995l) 



Delaunay Tessellation 

Construction of the Delaunay tessellation from the p oint sample (see fig. [TTI) . 
While we still use our own Voronoi-Delaunay code dvan de Wevgaerli 119941) . 
at present there is a score of efficient library routines available. Particularly 
noteworthy is the CGAL initiative, a large library of computational geometry 
routine^ 



''CGAL is a C++ library of algorithms and data structures for Computational Geometry, see 
[www. cgal.org I 
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• Field values point sample 

Dependent on whether it concerns the densities at the sample points or a mea- 
sured field value there are two options: 
+ General (non-density) field: 

(Sampled) value of field at sample point. 

+ Density field: The density values at the sampled points are determined from 
the corresponding Voronoi tessellations. The estimate of the density at each 
sample point is the normalized inverse of the volume of its contiguous 
Voronoi cell "W, of each point ;. The contiguous Voronoi cell of a point 
i is the union of all Delaunay tetrahedra of which point / forms one of the 
four vertices (see fig.|28]for an illustration). We recognize two applicable 
situations: 



- uniform sampling process: the point sample is an unbiased sample of 
the underlying density field. Typical example is that of A^-body simulation 
particles. For D-dimensional space the density estimate is, 

with Wi the weight of sample point /, usually we assume the same "mass" 
for each point. 

- systematic non-uniform sampling process: sampling density according to 
specified selection process quantified by an a priori known selection func- 
tion tl/{x), varying as function of sky position {a, 5) as well as depth/redshift. 
For D-dimensional space the density estimate is , 

p(x,) = (1-hD)— -. (17) 



Field Gradient 

Calculation of the field gradient estimate V/|„, in each D-dimensional Delau- 
nay simplex m(D = 3: tetrahedron; D -2: triangle) by solving the set of linear 
equations for the field values at the positions of the (Z) + 1) tetrahedron vertices, 



V/|„ 



/() /i fi h 



ro ri ¥2 r3 



(18) 



Evidently, linear interpolation for a field / is only meaningful when the field 
does not fluctuate strongly. Particularly relevant for velocity field reconstruc- 
tions is that there should be no orbit crossing flows within the volume of the 
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Delaunay cell which would involve multiple velocity values at any one loca- 
tion. In other words, DTFE velocity field analysis is only significant for laminar 
flows. 

Note that in the case of the sampled field being the velocity field v we may not 
only infer the velocity gradient in a Delaunay tetrahedron, but also the directly 
related quantities such as the velocity divergence, shear and vorticity. 



VTFE: 



DTFE: 



grid-based: 




Fig. 24. Two-dimensional display-grids in the VTFE, DTFE and grid-based reconstruction 
methods. The grid is overlaid on top of the basic underlying structure used to reconstruct 
the density field. SPH-like methods are not shown, because of the inherent difficulty in vi- 
sualizing their underlying structure, which does not consist of a subdivision of space in dis- 
tinct non-overlapping structural elements, but of circles of different radius at each position in 
space. 



• Interpolation. 

The final basic step of the DTFE procedure is the field interpolation. The pro- 
cessing and postprocessing steps involve numerous interpolation calculations, 
for each of the involved locations x. Given a location x, the Delaunay tetrahe- 
dron m in which it is embedded is determined. On the basis of the field gradient 
V/lm the field value is computed by (linear) interpolation (see eq. [TsT l. 

/(x) = /(X,) + V/[^_-(x-x,). (19) 

In principle, higher-order interpolation procedures are also possible. Two rele- 
vant high-order procedures are; 

- Spline Interpolation 

- Natural Neighbour Interpolation (see eqn.|8]and eqn.|9]l. 
Implementation of natural neighbour interpolation, on the basis of CGAL rou- 
tines, is presently in progress. 

• Processing. 

Though basically of the same character for practical purposes we make a dis- 
tinction between straightforward processing steps concerning the production of 
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Fig. 25. Image of a characteristic filamentary region in the outcome of a cosmological A'^- 
body simulation. Left-hand frame: particle distribution in a thin slice through the simulation 
box. Right-hand frame: tw o-dimensional slice through the three-dimensional DTFE density 
field reconstruction. From i lSchaatil2007h . 

images and simple smoothing filtering operations on the one hand, and more 
complex postprocessing on the other hand. The latter are treated in the next 
item. Basic to the processing steps is the determination of field values follow- 
ing the interpolation procedure(s) outlined above. 

Straightforward "first line" field operations are "Image reconstruction" and, 
subsequently, "Smoothing/Filtering ". 

+ Image reconstruction. 

For a set of image points, usually grid points, determine the image value: 
formally the average field value within the corresponding gridcell. In prac- 
tice a few different strategies may be followed, dictated by accuracy re- 
quirements. These are: 

- Formal geometric approach: integrate over the field values within each 
gridcell. This implies the calculation of the intersection of the relevant 
Delaunay tetrahedra and integration of the (linearly) running field values 
within the intersectio. Subsequently the integrands of each Delaunay inter- 
section are added and averaged over the gridcell volume. 

- Monte Carlo approach: approximate the integral by taking the average 
over a number of (interpolated) field values probed at randomly distributed 
locations within the gridcell around an image point. Finally average over 
the obtained field values within a gridcell. 

- Singular interpolation approach: a reasonable and usually satisfactory al- 
ternative to the formal geometric or Monte Carlo approach is the shortcut 
to hmit the field value calculation to that at the (grid) location of the im- 
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age point. This offers a reasonable approximation for gridcells which are 
smaller or comparable to that of intersecting Delaunay cells, on the condi- 
tion the field gradient within the cell(s) is not too large. 
+ Smoothing and Filtering: 

- Linear filtering of the field /: convolution of the field / with a filter func- 
tion Ws{x, y), usually user-specified. 



/,(x) 



f(x')W,{x',y)dx' 



(20) 



- Median (natural neighbour) filtering: the DTFE density field is adaptively 
smoothed on the bais of the median value of densities within the contiguous 
Vor onoi cell, the region defined by the a po int and its natural neighbours 
(see lPlaten. van de Weygaert & Jonesll2007h . 

- (Nonlinear) diffiision filtering: filtering of (sharply defined) fea tures in 



image s by solving appropriately defined diffusion equations (see e.g lMitra & Sicuranza , 
200d) . 








a 


■ 





Fig. 26. Filtering of density fields. The left-hand frame depicts the original DTFE density 
field. The subsequent frames show filtered DTFE density fields. The FWHM of the Gaussian 
filter is indicated by the shaded circle in the lower left-hand comer of these frames. From 
iSchaap, 2007). 



• Post-processing. 

The real potential of DTFE fields may be found in sophisticated applications, 
tuned towards uncovering characteristics of the reconstructed fields. An impor- 
tant aspect of this involves the analysis of structures in the density field. Some 
notable examples are: 



+ Advanced filtering operations . Potentially interestin g applications are those 



based on the use of wavelets dMartmez et al.L l2005h 



Cluster, Fi lament and Wall detection by means of the Multiscale Morphol- 
ogy Filter dAragon-Calvol 120071: lAragon-Calvo et alll2007l) . 



Void identification on the basis of the cosmic watershed algorithm dPlaten. van de Weygaert & Jones , 



2003). 



+ Halo detection in N-body simulations dNevrinck. Gnedin & Hamiltonll2005h . 
-I- The computation of 2-D surface densities for the study of gravitational lens- 



ing dBradac et al.L 120041) 
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In addition, DTFE enables the simultaneous and combined analysis of density 
fields and other relevant physical fields. As it allows the simultaneous determi- 
nation of density and velocity fields, it can serve as the basis for studies of the 
dynamics of structure formation in the cosmos. Its ability to detect substructure 
as well as reproduce the morphology of cosmic features and objects implies 
DTFE to be suited for assessing their dynamics without having to invoke artifi- 
cial filters. 



-I- DTFE as basis for the study of the full phase-space structure of structures 
and objects. The phase-space structure dark haloes in cosmolo gical struc- 
ture formation scenarios has been studied bv lArad et al.l (120041) . 



10 DTFE: densities and velocities 




Fig. 27. Illustration of a contiguous Voronoi cell. A contiguous Voronoi cell is the union of 
all Delaunay tetrahedra of which a point / is one of the vertices. 



10.1 DTFE density estimates 

The DTFE procedure extends the concept of interpolation of field values sampled 
at the point sample !P to the estimate of the density p(x) from the spatial point 
distribution itself. This is only feasible if the spatial distribution of the discrete point 
sample forms a fair and unbiased reflection of the underlying density field. 

It is commonly known that an optimal estimate for the spatial density at the loca- 
tion of a point x, in a discrete point sam ple P is given by th e inverse of the volume 
of the corresponding Voronoi cell (see Okabe et al. , 2000i section 8.5, for refer- 
ences). Tessellation-based methods for estimating the density have been introduced 
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by iBrownl ( Il965h and lOrdI ( 119781) . In astronomy, lEbeling & WiedenmannI ( Il993h 
were the first to use tessellation-based density estimators for the specific purpose of 
devising source detection al gorithms. This work has recently been applied t o clus 



ter detection algorithms b y (IRamella et al.L 12001 



2002t iLopes et al.L 120041) . Along the same lines, 



Kim et al.L 12002; Marino ni et al 



Ascasibar & BinnevI ( 12005 ) sug- 



gested that the use of a multidimensional binary tree might offer a computationally 
more efficient alternative. However, these studies have been restricted to raw esti- 
mates of the local sampling density at the position of the sampling points and have 
not yet included the more elaborate interpolation machinery of the DTFE and Nat- 
ural Neighbour Interpolation methods. 



Density definition 

The density field reconstruction of the DTFE procedure consists of two steps, the 
zeroth-order estimate po of the density values at the location of the points in P and 
the subsequent linear interpolation of those zeroth-order density estimates over the 
corresponding Delaunay grid throughout the sample volume. This yields the DTFE 
density field estimate p(x). 

It is straightforward to infer (see next sect. 110.11 ) that if the the zeroth-order es- 
timate of the density values would be the inverse of the regular Voronoi volume the 
condition of mass conservation would not be met. Instead, the DTFE procedure em- 
ploys a slightly modified yet related zeroth-order density estimate, the normalized 
inverse of the volume Vi^i) of the contiguous Voronoi cell 'W, of each point /. For 
Z)-dimensional space this is 

p(x,) = (1 +D) '—. (21) 

The contiguous Voronoi cell of a point ; is the union of all Delaunay tetrahedra of 
which point / forms one of the four vertices (see fig. |28] l. 



Mass Conservation 

An essential requirement for the cosmological purposes of our interpolation scheme 
is that the estimated DTFE density field p(x) should guarantee mass conservation: 
the total mass corresponding to the density field should be equal to the mass repre- 
sented by the sample points. Indeed, this is an absolutely crucial condition for many 
applications of a physical nature. Since the mass M is given by the integral of the 
density field p(x) over space, this translates into the integral requirement 




N 

= Y^nii = M = est., (22) 

1=1 
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Fig. 28. Relation between density and volume contiguous Voronoi cells. Two example points 
embedded within a filamentary structure (see fig.ll7t. 

with m, = m is the mass per sample point: the interpolation procedure should con- 
serve the mass M. 

The integral eq.|22]is equal to the volume below the linearly varying p-surface 
in (x,p)-space. In this space each Delaunay tetrahedron m is the base "hyper-plane" 
of a polyhedron 23,^, (here we use the term "tetrahedron" for any multidimensional 
Delaunay simplex). The total mass corresponding to the density field may therefore 
be written as the sum of the volumes of these polyhedra. 
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Nt 



M = J]V(&J, (23) 

m=l 

with Nt being the total number of Delaunay tetrahedra in the tessellation and V{DU 
the volume of polyhedron . This volume may be written as the average density 
at the vertices of Delaunay tetrahedron £)„, times its volume V(Dm) 

V(D*J = (a«i +Pna + ■■■ +P„,(Z)+i)) ViD,„) . (24) 

The points [ml, ml, . . ., m{D + 1)) are the nuclei which are vertices of the Delaunay 
tetrahedron 2)*,. The total mass M contained in the density field is the sum over all 
Delaunay tetrahedra within the sample volume: 

J Nt 

M = -^—^ ^ (p„,i +p,„2 + . . . +P,„(D+1)) V{Dm) . (25) 

A simple reordering of this sum yields 



in which D,„ i is one of the Noj Delaunay tetrahedra of which nuclues ; is a vertex. 
The complete set 2)„,_; constitutes teh contiguous Voronoi cell 'W, of nucleus /. Mass 
M may therefore be written as 

1 

M = ^—^ ^p, y(ni^,). (27) 

/=i 

This equation immediately implies that M is equal to the total mass of the sampling 
points (eq.l22li on the condition that the density at the location of the sampling points 
is 

p(x,) = (D + 1) '— . (28) 

This shows that the DTFE density estimate (eq. 1211 ). proportional to the inverse of 
contiguous Voronoi cell volume, indeed guarantees mass conservation. The corre- 
sponding normalization constant (1 + D) stands for the number of times each De- 
launay tetrahedron is used in defining the DTFE density field, equal to the number 
of ( 1 +D) sample points constituting its vertices. 



Non-uniform sampling 

Astronomical situations often involve a non-uniform sampling process. Often the 
non-uniform selection may be quantified by an a priori known selection function 
^/'(x). This situation is typical for galaxy surveys: i/r(x) may encapsulate differences 
in sampling density as function of 
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• Sky position {a, 6) 

In practice, galaxy (redshift) surveys are hardly ever perfectly uniform. Various 
kinds of factors - survey definition, observing conditions, instrumental aspects, 
incidental facts - will result in a non-uniform coverage of the sky. These may be 
encapsulated in a survey mask (//{a, 5). 

• Distance/redshift: (/';('") 

Magnitude- or flux-limited surveys induce a systematic redshift selection. At 
higher redshifts, magnitude-limited surveys can only probe galaxies whose lu- 
minosity exceeds an ever increasing value L{z). The corresponding radial selec- 
tion function ip- is given by 

Uz) = , (29) 

L '^^^^^^ 

where 0{L) is the galaxy luminosity function. 

Both selection effects will lead to an underestimate of density value. To correct for 
these variations in sky completeness and/or redshift/depth completeness the esti- 
mated density at a sample point/galaxy ; is weighted by the inverse of the selection 
function, 

ilf{iLi) ^ ilfs{ai,6i)il,,{ri), (30) 

so that we obtain, 

in ' 

. While it is straightforward to correct the density estimates accordingly, a com- 
pUcation occurs for locations where the completeness is very low or even equal to 
zero. In general, regions with redshift completeness zero should be excluded from 
the correction procedure, even though for specific cosmological contexts it is fea- 
sible to incorporate field reconstruc tion procedures utihzing field correlations. A 
constrained random field approach ( Bertschinger . 1987 : Hofiman & RibakL 1991 



Zaroubi et al.Lll995h uses the autocorrelation function of the presum ed density field 



to inf er a realization of the corresponding Gaussian field. Recently, ( Aragon-Calvoi 



20071) developed a DTFE based technique which manages to reconstruct a structured 
density field pattern within selection gaps whose size does not exceed the density 
field's coherence scale. 



10.2 DTFE Velocity Fields 

DTFE density and velocity field gradients 

The value of the density and velocity field gradient in each Delaunay tetrahedron is 
directly and uniquely determined from the location r = (x, y, z) of the four points 
forming the Delaunay tetrahedra's vertices, ro, ri, r2 and r3, and the value of the 
estimated density and sampled velocities at each of these locations, (po, vo), (p\,\i), 
(P2,V2) and(p3,V3), 



62 



Rien van de Weygaert & Willem Schaap 



PO Pi P2 P3 
Vo Vi V2 V3 

ro ri r2 r3 



(32) 



The four vertices of the Delaunay tetrahedron are both necessary and sufficient for 
computing the entire 3x3 velocity gradient tensor dvi/dxj. Evidently, the same 
holds for the density gradient 5p/5xy. We define the matrix A is defined on the basis 
of the vertex distances {Ax„,Ay„,Azn) (n=l,2,3), 



AXn 



X„ - Xq 

yn - yo 

Zn - ZO 



Axi Ayi Azi 
Ax2 Ay2 Az2 
Axi Ayi Azj. 



(33) 



Similarly defining z)v„ = v„-Vo(n = 1, 2, 3) and zip,, = p„ - po(n - 1,2,3) it 
is straightforward to compute directly and simultaneously the density field gradient 
Vp|„, and the velocity field gradient Vv|„, - dvi/dxj in Delaunay tetrahedron m via 
the inversion. 



dp \ 
dx 

dp 

dy 

dp 
dz 



= A 



Api 

Ap2 

Api. 



(dv. 


dVy 




dx 


dx 


dx 




dvy 


dv- 


dy 


'dy 


dy 


5vv 


dVy 


dv. 


dz 


"3? 


dz> 



= A 



(A\'x, Aviy A\'i-\ 

Av2x Av2y Av2z 
yV3.v ZfV3v Av^-J 



34) 



Velocity divergence, shear and vorticity 

From the nine velocity gradient components dvi/dxj we can directly determine the 
three velocity deformation modes, the velocity divergence V ■ v, the shear cr,j and 
the vorticity w, 

„ [dv, dvy dv, 
V-v - I — + — + ^ 

dx dy dz 



1 j dvi dvj 

2 [dxj dxi 

1 f dvi dvj 



3(V-v)5,v, 



(35) 



2 \dxj dxi 



where w = V x v = e^cjij (and is the completely antisymmetric tensor). In the 
theory of gravitational instability, there will be no vorticity contribution as long as 
there has not been shell crossing (ie. in the linear and mildly nonlinear regime). 
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10.3 DTFE: velocity field Umitations 

DTFE interpolation of velocity field values is only feasible in regions devoid of 
multistream flows. As soon as there are multiple flows - notably in high-density 
cluster concentrations or in the highest density realms of the filamentary and planar 
caustics in the cosmic web - the method breaks down and cannot be apllied. 

In the study presented here this is particularly so in high-density clusters. The 
complication can be circumvented by filtering the velocities over a sufficiently large 
region, imposing an additional resolution constraint on the DTFE velocity field. 
Implicitly this has actually already been accomplished in the linearization procedure 
of the velocity fields preceding the DTFE processing. The linearization of the input 
velocities involves a kernel size of Vs/i 'Mpc, so that the resolution of the velocity 
field is set to a lower limit of Vs/i 'Mpc. This is sufficient to assure the viability of 
the DTFE velocity field reconstructions. 

10.4 DTFE: mass-weighted vs. volume-weighted 

A major and essential aspect of DTFE field estimates is that it concerns volume- 
weighted field averages 



where W(x, y) is the adopted filter function defining the weight of a mass element 
in a way that is dependent on its position y with respect to the position x. 

Analytical calculations of physical systems in an advanced stage of evolution 
do quite often concern a perturbation analysis. In a cosmological context we may 
think of the nonlinear evolution of density and velocity field perturbations. In or- 
der to avoid unnecessary mathematical complications most results concern volume- 
weighted quantities. However, when seeking to relate these to the results of obser- 
vational results or numerical calculations, involving a discrete sample of measure- 
ment points, nearly all conventional techniques implicitly involve mass-weighted 
averages. This may lead to considerable confusion, even to wrongly motivated con- 
clusions. 

Conventional schemes for translating a discrete sample of field values fi into 
a continuous field /(x) are usually based upon a suitably weighted sum over the 
discretely sampled field values, involving the kernel weight functions W(x, y). It is 




(36) 



instead of the more common mass-weighted field averages. 




(37) 
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Straightforward to convert this discrete sum into an integral over Dirac delta func- 
tions, 



/(X) = 



X.^^/;mx-x,) 

J Jy^x-y) Y^.^^ My-x,) 

/ 



(38) 



dypiy) W(x-y) 



Evidently, a weighted sum implicitly yields a mass-weighted field average. Notice 
that this is true for rigid grid-based averages, but also for spatially adaptive SPH like 
evaluations. 

A reasonable approximation of volume-averaged quantities may be obtained by 
volume averaging over quantities that were mass-filtered with, in c omparison, a very 



small scale for the mass-weighting filter function. This prodded lJuszkiewicz et al 



(119951) to suggest a (partial) remedy in the form of a two-step scheme for velocity 
field analysis. First, the field values are interpolated onto a grid according to eqn.[39l 
Subsequently, the resulting grid of field values is used to determine volume-averages 
according to eqn.[36] We can then make the interesting observation that the asymp- 
totic limit of this, using a filter with an infinitely small filter length, yields 

N 



VwJixd /(Xl) + y ^/(X,) 



/„,,,(xo) = ^.^^ = ^ /(Xl) (39) 



,=2 ^1 

where we have ordered the locations ; by increasing distance to Xq and thus by 
decreasing value of w, . 

The interesting conclusion is that the resulting estimate of the volume-weighted 
average is in fact the field value at the location of the closest sample point Xi, /(xi). 
This means we should divide up space into regions consisting of that part of space 
closer to a particular sample point than to any of the other sample points and take 
the field value of that sample point as the value of the field in that region. This is 
nothing but dividing up space according to the Voronoi tessellation of the set of 
sample points P. This observation formed the rationale behind the introduction and 
definitio n of Voronoi and Delaunay tessellation interpolation methods for velocity 
fields by lBernardeau & van de Weygaert (Il996h . 



Geometric Analysis Cosmic Web 



65 



While the resulting Voronoi Tessellation Field Estimator would involve a dis- 
continuous field, the step towards a multidimensional linear interpolation procedure 
would guarantee the continuity of field values. The subsequent logical step, invok- 
ing the dual Delaunay tessellation as equally adaptive spatial linear interpolation 
grid, leads to the definition of the DTFE interpolation scheme. 



10.5 DTFE Density and Velocity Maps: non-uniform resolution 

For a proper quantitative analysis and assessment of DTFE density and velocity 
field reconstructions it is of the utmost importance to take into account the inhomo- 
geneous spatial resolution of raw DTFE maps. 

Cosmic density and velocity fields, as well as possible other fields of relevance, 
are composed of contributions from a range of scales. By implicitly filtering out 
small-scale contributions in regions with a lower sampling density DTFE will in- 
clude a smaller range of spatial scales contributing to a field reconstruction. Even 
while selection function corrections are taken into account, the DTFE density and 
velocity fields will therefore involve l ower amplitudes. DTFE velocity fields are 



less affected than DTFE density fields ( IRomano-Diaz & van de Wevsaerll l2007h . a 



manifestation of the fact that the cosmic velocity field is dominated by larger scale 
modes than the density field. 

In addition, it will also lead to a diminished "morphological resolution". The 
sampling density may decrease to such a level lower than required to resolve the ge- 
ometry or morphology of weblike features. Once this level has been reached DTFE 
will no longer be suited for an analysis of the Cosmic Web. 



11 DTFE: Technical Issues 

11.1 DTFE Kernel and Effective Resolution 

DTFE distributes the mass m, corresponding to each sampling point ; over its cor- 
responding contiguous Voronoi cell. A particular sampling point ; will therefore 
contribute solely to the mass of those Delaunay simplices of which it is a vertex. 
This restricted spatial smoothing is the essence of the strict locality of the DTFE 
reconstruction scheme. 

The expression for the interpolation kernel of DTFE provides substantial insight 
into its local spatial resolution. Here we concentrate on the density field kernel, the 
shape and scale of the interpolation kernels of other fields is comparable. 

In essence, a density reconstruction scheme distributes the mass m, of a sampling 
particle / over a region of space according to a distributing function Tiix), 



N 
;=1 



(40) 
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Fig. 29. Typical interpolation kernels for points embedded within three different morpholog- 
ical environments: a point in a high-density peak (bottom), a point in a filament (centre) and 
a point in a low-density field environment (top). The three kernels depicted concern three 
different environments: a point embedded within a Gaussian density peak (lefthand), a point 
within an elongated filamentary concentration of points (centre) and a point in a low-density 
environment (righthand). 
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with J dx Tiix) = 1 V /. Generically, the shape and scale of the interpolation kernel 
Tiix) may be different for each sample point. The four examples of nn-neighbour 
kernels in fig.|2T|are a clear illustration of this. 

For the linear interpolation of DTFE we find that (see eqn. fTTT i for a position x 
within a Delaunay tetrahedron m defined by {l+D) sample points {x,„o,x,„i, . . . ,x.„,o] 
the interpolation indices (f)dt,i are given by 



/(X) = 



\ + ti + t2 + . . . + to / e {mO, ml, ... , mD] 

(41) 

/ i {mO, ml, . . ., mD} 



In this, for / e {mO, ml, . . ., mD}, the D parameters {ti, . . ., to} are to be solved from 

X = X; + ri(x,„i - X,) + f2(x,„2 - X,) + . . . + toiXmD - X/) , (42) 

On the basis of eqns. ( |40] | and ( fTTT ) the expression for the DTFE kernel is easily 
derived: 



Yf^'PdtA^) if xe'W; 

(43) 

if xi^Vi 



in which TV,- is the contiguous Voronoi cell of sampling point /. In this respect we 
should realize that in two dimensions the contiguous Voronoi cell is defined by on 
average 7 sample points: the point itself and its natural neighbours. In three di- 
mensions this number is, on average, 13.04. Even with respect to spatially adaptive 
smoothing such as embodied by the kernels used in Smooth Particle Hydrodynam- 
ics, defined a certain number of nearest neighbours (usually in the order of 40), the 
DTFE kernel is indeed highly localized. 

A set of DTFE smoothing kernels is depicted in fig.|29] Their local and adaptive 
nature may be best appreciated from the comparison with the corresponding kernels 
for regular (rigid) gridbased TSC interpolation, a scale adaptive SPH smoothing 
kernel (based on 40 nearest neighbours) and a zeroth-order Voronoi (VTFE) ker- 
nel (where the density at a point is set equal to the inverse of the volume of the 
corresponding Voronoi cell). The three kernels depicted concern three different en- 
vironments: a point embedded within a Gaussian density peak (lefthand), a point 
within an elongated filamentary concentration of points (centre) and a point in a 
low-density environment (righthand). The figure clearly illustrates the adaptivity in 
scale and geometry of the DTFE kernel. 



11.2 Noise and Sampling characteristics 

In order to appreciate the workings of DTFE one needs to take into account the 
noise characteristics of the method. Tessellation based reconstructions do start from 
a discrete random sampling of a field. This will induce noise and errors. Discrete 
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Fig. 30. Poisson sampling noise in uniform fields. The rows illustrate three Poisson point 
samplings of a uniform field with increasing sampling density (from top to bottom consisting 
of 100, 250 and 1000 points). From left to right the point distribution, the corresponding 
Voronoi tessellation, the zeroth-order VTFE reconstructed density field, the first-order DTFE 
reconstructed density field and the corresponding Delaunay tessellation are shown. From 
Schaap 2007. 



sampling noise will propagate even more strongly into DTFE density reconstruc- 
tions as the discrete point sample itself will be the source for the density value 
estimates. 

In order to understand how sampling noise affects the reconstructed DTFE fields 
it is imperative to see how intrinsically uniform fields are affected by the discrete 
sampling. Even though a uniform Poisson process represents a discrete reflection of 
a uniform density field, the stochastic nature of the Poisson sample will induce a 
non-uniform distribution of Voronoi and Delaunay volumes (see sect. 16.5b . Because 
the volumes of Voronoi and/or Delaunay cells are at the basis of the DTFE density 
estimates their non-uniform distribution propagates immediately into fluctuations in 
the reconstructed density field. 

This is illustrated in fig. 111.11 in which three uniform Poisson point samples, 
each of a different sampling density «, are shown together with the correspond- 
ing Voronoi and Delaunay tessellations. In addition, we have shown the first-order 
DTFE reconstructed density fields, along with the zeroth-order VTFE density field 
(where the density at a point is set equal to the inverse of the volume of the cor- 
responding Voronoi cell). The variation in both the Delaunay and Voronoi cells, as 
well as in the implied VTFE and DTFE density field reconstructions, provide a rea- 
sonable impression of the fluctuation s going along with these interpolation schemes. 



Following the Kiang suggestion (lKiangill966l) for the Gamma type volume dis 



tribution for the volumes of Voronoi cells (eqn. lU, we may find an analytical ex- 
pression for the density distribution for the zeroth-order VTFE field; 



(44) 

p-«exp(-|) 3D 



1944 
5 



The normalized density value p is the density estimate p in units of the average 
density {p). While in principle this only concerns the zeroth-order density estimate, 
it turns out that these expression also provide an excellent description of the one- 
point distribution function of the first-order DTFE density field, both in two and 
three dimensions. This may be appreciated from Fig. [31] showing the pdfs for DTFE 
density field realizations for a Poisson random field of 10,000 points (2-D) and 
100,000 point (3-D). The superimposed analytical expressions (eqn. l44li represent 
excellent fits. 

The 2-D and 3-D distributions are clearly non-Gaussian, involving a tail extend- 
ing towards high density values. The positive value of the skewness (2 V3 for 2D 




and Vs for 3D) confirms the presence of this tail. Interestingly, the distribution falls 
off considerably more rapid for d - 3 than for d - 2. Also we see that the distribu- 
tions are more narrow than in the case of a regular Gaussian: the variance is 1 /3 for 
d - 2 and 1/5 for d=3, confirmed byt the strongly positive value of the curtosis. The 
larger value for d - 2 indicates that it is more strongly peaked than the distribution 
forc/ = 3. 

On the basis of the above one may also evaluate the significance of DTFE den- 
sity field reconstructio ns, even includ ing that for intrinsically inhomogeneous fields. 
For details we refer to Schaapl ( 2007 ). 



11.3 Computational Cost 

. The computational cost of DTFE is not overriding. Most computational effort is 
directed towards the computation of the Delaunay tessellation of a point set of 
particles. The required CPU time is in the order of 0{N log N), comparable to the 
cost of generating the neighbour list in SPH. The different computational steps of 
the DTFE procedure, along with their scaling as a function of number of sample 
points A^, are: 

1 . Construction of the Delaunay tessellation: 0{N log A'); 

2. Construction of the adjacency matrix: 0(Ny, 

3. Calculation of the density values at each location: 0{Ny, 

4. Calculation of the density gradient inside each Delaunay triangle: 0(N); 

5. Interpolation of the density to an image grid: 0(ngrid^ ■ N^'^). 
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Step2, involving the calculation of the adjacency matrix necessary for the walking 
triangle algorithm used for Delaunay tetrahedron identification, may be incorpo- 
rated in the construction of the Delaunay tessellation itself and therefore omitted. 
The last step, the interpolation of the density to an image grid, is part of the post- 
processing operation and could be replaced by any other desired operation. It mainly 
depends on the number of gridcells per dimension. 



12 DTFE: Hierarchical and Anisotropic Patterns 

To demonstrate the ability of the Delaunay Tessellation Field Estimator to quanti- 
tatively trace key characteristics of the cosmic web we investigate in some detail 
two aspects, its ability to trace the hierarchical matter distribution and its ability to 
reproduce the shape of anisotropic - filamentary and planar - weblike patterns. 



12.1 Hierarchy and Dynamic Range 



The fractal Soneira-Peebles model ( ISoneira & Peeblesl 119771) has been used to as 



sess the level to which DTFE is able to trace a density field over the full range of 
scales represented in a point sample. The Soneira-Peebles model is an analytic self- 
similar spatial point distribution which was defined for the purpose of modelling the 
galaxy distribution, such that its statistical properties would be tuned to reality. An 
important property of the Soneira-Peebles model is that it is one of the few nonlin- 
ear models of the galaxy distribution whose statistical properties can be fully and 
analytically evaluated. This concerns its power-law two-point correlation function, 
correlation dimension and its Hausdorff dimension. The Soneira-Peebles model is 
specified by three parameters. The starting point of the model is a level-0 sphere 
of radius R. At each level-;w a number of rj subspheres are placed randomly within 
their parent level-m sphere: the level-(m +1) spheres have a radius R/A where A > I, 
the size ratio between parent sphere and subsphere. This process is repeated for L 
successive levels, yielding r]^ level-L spheres of radius R/A^. At the center of each 
of these spheres a point is placed, producing a point sample of 77^ points. While 
this produces a pure singular Soneira-Peebles model, usually a set of these is su- 
perimposed to produce a somewhat more realistically looking model of the galaxy 
distribution, an extended Soneira-Peebles model. 

An impression of the self-similar point sets produced by the Soneira-Peebles 
process may be obtained from fig. 112. II The top row contains a series of three point 
distributions, zoomed in at three consecutive levels on a singular Soneira-Peebles 
model realization with (77 = 4, /I = 1.90, L - 8). The next three columns depict 
the density field renderings produced by three different interpolation schemes, a 
regular rigid grid-based TSC scheme, a spatially adaptive SPH scheme and finally 
the DTFE reconstruction. The figure clearly shows the supreme resolution of the 
latter. By comparison, even the SPH reconstructions appear to fade for the highest 
resolution frame. 



Geometric Analysis Cosmic Web 



73 



Fig. 32. Singular Soneira-Peebles model with 77 = 4, /i = 1.9 and L = 8. Top row: full point 
distribution (left-hand frame) and zoom-ins focusing on a particular structure (central and 
right-hand frames). Rows 2 to 4: corresponding density field reconstructions produced using 
the TSC, SPH and DTFE methods. 



Self-simUarity 

One important manifestations of the self-similarity of the defined Soneira-Peebles 
distribution is reflected in the power-law two-point correlation function. For M di- 
mensions it is given by 

(45) 

y ^ M-l — ^ for — — <r<R. 



log A I A^-^ 

The parameters 77 and A may be adjusted such that they yield the desired value for 
the correlation slope y. 

To probe the self-similarity we look at the one-point distribution of the density 
field, both for the point distribution as well as the TSC, SPH and DTFE density 
field reconstructions. Mathematically, the condition of self-similarity implies that 
the PDF corresponding to a density field p(x) inside an n-level circle of radius R/A" 
should be equal to the PDF inside the reference circle of radius R, after the density 
field in the n-level circle has been scaled according to 

p(x) ^ p„(x) = p(x)M*^'', (46) 

in which M is the dimension of space. Yet another multiplication factor of A^" has to 
be included to properly normalize the PDF (per density unit). In total this results in 
an increase by a factor A-'^". Self-similar distributions would look exactly the same 
at different levels of resolution once scaled accordingly. This self-similarity finds its 
expression in a universal power-law behaviour of the PDFs at different levels. 

We have tested the self-similar scaling of the pdfs for a range of Soneira-Peebles 
models, each with a different fractal dimensions dSchaap & van de Wevgaert , 



For a Soneira-Peebles model with parameters {j] - 2,A - 1.75, L = 14) we show 
the resulting scaled PDFs for the TSC, SPH and DTFE density field reconstructions 
in the top row of Fig. [33]. The self-similar scaling of the TSC density field recon- 
structions is hardly convincing. On the other hand, while the SPH reconstruction 
performs considerably better, it is only the DTFE rendering which within the den- 
sity range probed by each level displays an almost perfect power-law scaling ! Also 
notice that the DTFE field manages to probe the density field over a considerably 
larger density range than e.g. the SPH density field. Subsequently we determined 
the slope a of the "universal " power-law PDF and compared it with the theoretical 
predictions. The set of three frames in the bottom row of Fig. [33] show the result- 
ing distributions with the fitted power-laws. Going from left to right, the frames 
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Fig. 33. Top row: Average PDFs of the density field in circles of different level (see text 
for a description) for the TSC, SPH and DTFE density field reconstructions. Model with 
Tj = 2, A = 1.75 and L = 14. Bottom row: scaled PDFs of Soneira-Peebles density field 
reconstructions. Each frame corresponds to a Soneira-Peebles realization of a difi"erent fractal 
dimension, denoted in the upper right-hand corner. In each frame the TSC, SPH and DTFE 
reconstructed PDFs are shown. 



D 


Q-(theory) 


a(TSC) 


a(SPU) 


Q-CDTFE) 


0.63 


-1.69 


-0.81 


-1.32 


-1.70 


0.86 


-1.57 


-0.82 


-1.24 


-1.60 


1.23 


-1.39 


-0.79 


-1.13 


-1.38 



Table 1. Slopes of the power-law region of the PDF of a Soneira-Peebles density field as 
reconstructed by the TSC, SPH and DTFE procedures. The theoretical value (Eqn.l48b is also 
listed. Values are listed for three different Soneira-Peebles realizations, each with a different 
fractal dimension D. 



in this figure correspond to Soneira-Peebles realizations with fractal dimensions of 
D = 0.63, D = 0.85 and D = 1.23. The slope a of the PDF can be found by 
comparing the PDF at two different levels, 

^ logPDF(pi) - logPDF(po) 
logpi - logpo 

(47) 

_ \og{A^'^"/j]") _ D 2 
" log(l/i'^'') ~ M ~ ' 

in which D is the fractal dimension of the singular Soneira-Peebles model. When 
turning to table [T] we not only find that the values of a derived from the TSC, 
SPH and DTFE fields do differ significantly, a fact which has been clear borne out 
by fig. |33] but also that the DTFE density field PDFs do reproduce to an impres- 
sively accurate level the analytically expected power-law slope of the model itself 



dSchaap & van de Wevgaerlll2007bl) . It is hardly possible to find a more convincing 



argument for the ability of DTFE to deal with hierarchical density distributions ! 



12.2 Shapes and Patterns 



DTFE's ability to trace anisotropic weblike patterns is tested on the basis of a 
class of heuristic models of cellular matter distributions, Voron oi Clustering Models 
dSchaap & van de Wevgaertl l2007d; Ivan de Wevgaern, l2007ah . Voronoi models use 
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the frame of a Voronoi tessellation as the skeleton of the galaxy distribution, iden- 
tifying the structural frame around which matter will gradually assemble during the 
emergence of cosmic structure. The interior of Voronoi cells correspond to voids and 
the Voronoi planes with sheets of galaxies. The edges delineating the rim of each 
wall are identified with the filaments in the galaxy distribution. What is usually de- 
noted as a flattened "supercluster" will comprise an assembly of various connecting 
walls in the Voronoi foam, as elongated "superclusters" of "filaments" will usually 
consist of a few coupled edges. The most outstanding structural elements are the 
vertices, corresponding to the very dense compact nodes within the cosmic web, the 
rich clusters of galaxies. 

The Voronoi clustering models offer flexible templates for cellular patterns, and 
they are easy to tune towards a particular spatial cellular morphology. To investigate 
the shape performance of DTFE we use pure Voronoi Element Models, tailor-made 
heuristic "galaxy" distributions either and exclusively in and around 1) the Voronoi 
walls, 2) the Voronoi edges, 3) the Voronoi vertices. Starting from a random initial 
distribution, all model galaxies are projected onto the relevant wall, edge or vertex 
of the Voronoi cell in which they are located. 

The three boxes in the top row of fig. [34l depict a realization for a 3-D Voronoi 
wall model, a Voronoi filament model and a Voronoi cluster model. Planar sections 
through the TSC, SPH and DTFE density field reconstructions of each of these mod- 
els are shown in three consecutive rows, by means of greyscale maps. The distinctly 
planar geometry of the Voronoi wall model and the one-dimensional geometry if 
the Voronoi filament model is clearly recognizable from the sharply defined features 
in the DTFE density field reconstruction. While the SPH reconstructions outline 
the same distinct patterns, in general the structural features have a more puffy ap- 
pearance. The gridbased TSC method is highly insensitive to the intricate weblike 
Voronoi features, often the effective rigid gridscale TSC filter is not able to render 
and detect them. 

The DTFE reconstruction of the Voronoi cluster models (fig. [34l lower right- 
hand) does depict some of the artefacts induced by the DTFE technique. DTFE 
is succesfull in identifying nearly every cluster, even the poor clusters which SPH 
cannot find. The compact dense cluster also present a challenge in that they reveal 
low-density triangular wings in the regions between the clusters. Even though these 
wings include only a minor fraction of the matter distribution, they are an arte- 
fact which should be accounted for. Evidently, the SPH reconstruction of individual 
clusters as spherical blobs is visually more appealing. 

Voronoi Filament Model 

The best testbench for the ability of the different methods to recover the patterns and 
morphology of the fully three-dimensional density field is that of the contrast rich 
Voronoi filament models. In Fig.[35]the central part of a sample box of the Vorotioi 
filament model realization is shown. Isodensity contour levels are chosen such that 
65% of the mass is enclosed within regions of density equal to or higher the corre- 
sponding density value. The resulting TSC, SPH and DTFE density fields are de- 
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Fig. 35. Three-dimensional visualization of the Voronoi filament model and the correspond- 
ing TSC, SPH and DTFE density field reconstructions. The density contours have been cho- 
sen such that 65% of the mass is enclosed. The an^ows indicate two structures which are 
visible in both the galaxy distribution and the DTFE reconstruction, but not in the TSC and 
SPH reconstructions. 



picted in the lower lefthand (TSC), lower righthand (SPH) and top frame (DTFE). 
The galaxy distribution in the upper lefthand frame does contain all galaxies within 
the region. Evidently, the galaxies have distributed themselves over a large range of 
densities and thus occupy a larger fraction of space than that outlined by the density 
contours. The appearances of the TSC, SPH and DTFE patterns do differ substan- 
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Fig. 36. Anisotropy measurements for the Voronoi models. Plotted is the longest-to-shortest 
axis ratio of the intertia tensor inside a series of concentric spheres centered on a characteristic 
structure as a function of the radius of the sphere. The radius is given in units of the standard 
deviation (cr) of the corresponding Gaussian density profiles. The lefthand frame corresponds 
to the Voronoi wall model, the central frame to the Voronoi filament model and the righthand 
frame to the Voronoi cluster model. In each frame the results are shown for the TSC, SPH and 
DTFE reconstructions, as well as for the galaxy distribution. The meaning of the symbols is 
depicted in the right-hand frame. 

tially. Part of this is due to a different effective scale of the filter kernel. The 65% 
mass contour corresponds to a density contour p - 0.55 in the TSC field, p = 1 .4 
in the SPH reconstruction and p = 2.0 in the DTFE reconstruction (p in units of the 
average density). The fine filamentary maze seen in the galaxy distribution is hardly 
reflected in the TSC grid based reconstruction even though the global structure, the 
almost ringlike arrangement of filaments, is clearly recognizable. The SPH density 
field fares considerably better, as it outlines the basic configuration of the filamen- 
tary web. Nonetheless, the SPH rendered filaments have a width determined by the 
SPH kernel scale, resulting in a pattern of tubes. Bridging substantial density gra- 
dients is problematic for SPH reconstructions, partly due to the lack of directional 
information. 

It is the DTFE reconstruction (top frame fig.[35]l which yields the most outstand- 
ing reproduction of the filamentary weblike character of the galaxy distribution. A 
detailed comparison between the galaxy distribution and the density surfaces show 
that it manages to trace the most minute details in the cosmic web. Note that the den- 
sity contours do enclose only 65% of the mass, and thus relates to a smaller volume 
than suggested by the features in the galaxy distribution itself. The success of the 
DTFE method is underlined by identifying a few features in the galaxy distribution 
which were identified by DTFE but not by SPH and TSC. The arrows in Fig. |35] 
point at two tenuous filamentary features visible in the galaxy distribution as well as 
in the DTFE field, yet entirely absent from the TSC and SPH fields. In comparison 
to the inflated contours of the SPH and TSC reconstructions, the structure outlined 
by the DTFE density field has a more intricate, even somewhat tenuous, appearance 
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marked by a substantial richness in structural detail and contrast. Some artefacts 
of the DTFE method are also visible: in particular near intersections of filaments 
we tend to find triangular features which can not be identified with similar struc- 
tures in the galaxy distribution. Nearby filaments are connected by relatively small 
tetrahedra, translating into high density features of such shape. 

Shape and Morphology Analysis 

An important measure of the local density distribution concerns the shape of the 
density contour levels. Various representative features in the three Voronoi element 
models were identified, followed by a study of their shape over a range of spa- 
tial scales. The axis ratios of the local density distribution, within a radius R, was 
computed from the eigenvalues of the mass intertia tensor The results of the shape 
analysis are shown in Fig. [36] From left to right, the three frames present the axis 
ratio of the longest over the smallest axis, a 1/03, for walls, filaments and clusters, 
as a function of the scale R over which the shape was measured. The open circles 
represent the shape of the particle distribution, the triangles the shape found in the 
equivalent DTFE density field, while crosses and squares stand for the findings of 
SPH and TSC. 

In the case of the Voronoi cluster models all three density field reconstructions 
agree with the sphericity of the particle distributions. In the central and righthand 
frame of Fig. [36] we see that the intrinsic shapes of the walls and filaments become 
more pronounced as the radius R increases. The uniform increase of the axis ratio 
a 1/(33 with R is a reflection of the influence of the intrinsic width of the walls and 
filaments on the measured shape. For small radii the mass distribution around the 
center of one of these features is largely confined to the interior of the wall or fila- 
ment and thus near-isotropic. As the radius R increases in value, the intrinsic shape 
of these features comes to the fore, resulting in a revealing function of shape as 
function of R. 

The findings of our analysis are remarkably strong and unequivocal: over the 
complete range of radii we find a striking agreement between DTFE and the cor- 
responding particle distribution. SPH reveals systematic and substantial differences 
in that they are tend to be more spherical than the particle distribution, in particu- 
lar for the strongly anisotropic distributions of the walls and filaments. In turn, the 
SPH shapes are substantially better than those obtained from the TSC reconstruc- 
tions. The rigidity of the gridbased TSC density field reconstructions renders them 
the worst descriptions of the anisotropy of the local matter distribution. These re- 
sults show that DTFE is indeed capable of an impressively accurate description of 
the shape of walls, filaments and clusters. It provides a striking confirmation of the 
earlier discussed visual impressions. 
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13 DTFE: Velocity Field Analysis 



De facto the DTFE method has been first defined in the context of a description and 
analysis of cosmic flow fields which have b een sampled by a set of discretely and 
sparsely sampled galaxy peculiar velocities. iBernardeau & van de Wevgaerll ( ll996h 
demonstrated the method's superior performance with respect to conventional inter- 
polation procedures. In particular, they also proved that the obtained field estimates 
involved those of the proper volume-weighted quantities, instead of the usually im- 
plicit mass-weighted quantities (see sect. 110.4b . This corrected a few fundamental 
biases in estimates of higher order velocity field moments. 





Fig. 37. The density and velocity field of the LCDM GIF N-body simulation, computed by 
DTFE. The bottom two frames show the density (bottom left) and velocity field (bottom 
right) in a central slice through the complete simulation box. The density map is a grayscale 
map. The velocity field is shown by means of velocity vectors: the vectors are the velocities 
component within the slice, their size proportional to the velocity's amplitude. By means of 
DTFE the images in both top frames zoom in on the density structure (left) and flow field 
(bottom) around a filamentary structure (whose location in the box is marked by the solid 
square in the bottom righthand panel). The shear flow along the filaments is meticulously 
resolved. Note that DTFE does not need extra operations to zoom in, one DTFE calculation 
suffices to extract all necessary information. From Schaap 2007. 
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The potential of the DTFE formahsm may be fully exploited within the context 
of the analysis of A^-body simulations and the galaxy distribution in galaxy redshift 
surveys. Not only it allows a study of the patterns in the nonlinear matter distribution 
but also a study of the related velocity flows. Because DTFE manages to follow both 
the density distribution and the corresponding velocity distribution into nonlinear 
features it opens up the window towards a study of the dynamics of the formation of 
the cosmic web and its corresponding elements. Evidently, such an analysis of the 
dynamics is limited to regions and scales without multistream flows (see sect. 110. 31 ). 

A study of a set of GIF ACDM simulations has opened up a detailed view of 
the dynamics in and around filaments and other components of the cosmic web by 
allowing an assessment of the detailed density and velocity field in and around them 
(see fig.l37Ti. DTFE density and velocity fields may be depicted at any arbitrary reso- 
lution without involving any extra calculation: zoom-ins represent themselves a real 
magnification of the reconstructed fields. This is in stark contrast to conventional 
reconstructions in which the resolution is arbitrarily set by the users and whose 
properties are dependent on the adopted resolution. 



13.1 A case study: Void Outflow 

In Fig. [38] a typical void-like region is shown, together with the resulting DTFE 
density and velocity field reco nstructions. It concerns a voidlike region selected 
from a ACDM GIF simulation ( IKaufiinann et al. . Il999h . Figure [To] shows a major 



part of the (DTFE) density field of the entire GIF simulation which contains the 
void. It concerns a 256^ particles GIF A^-body simulation, encompassing a ACDM 
(Qm = 0.3, Qa - 0.7, Ho = 70km/s/Mpc) density field within a (periodic) cubic 
box with length 141/7 'Mpc and produced by means of an adaptive P^M A^-body 
code. 

The top left frame shows the particle distribution in and around the void within 
this 42.5/2 'Mpc wide and l/j 'Mpc thick slice through the simulation box. The 
corresponding density field (top righthand frame) and velocity vector field (bottom 
lefthand frame) are a nice illustration of the ability of DTFE to translate the inho- 
mogeneous particle distribution into highly resolved continuous and volume-filling 
fields. 

The DTFE procedure clearly manages to render the void as a slowly varying 
region of low density. Notice the clear distinction between the empty (dark) interior 
regions of the void and its edges. In the interior of the void several smaller sub- 
voids can be distinguished, with boundaries consisting of low density filamentary or 
planar structures. Such a hierarchy of voids, with large voids containing the traces 
of the smaller ones from which it f ormed earlier through merging, has been de- 
scribed by theories of void evolution jRegos & Ge"iie3.ll99ll:lDubinski et 19931 



van de Weygaert & van KampenL 1993^ Sheth & van de Weygaertl 20041) 



The velocity vector field in and around the void represents a nice illustration of 
the qualities of DTFE to recover the general velocity flow. It also demonstrates its 
ability identify detailed features within the velocity field. The flow in and around 
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Fig. 38. The density and velocity field around a void in the GIF LCDM simulation. The top 
righthand panel shows the N-body simulation particle distribution within a slice through the 
simulation box, centered on the void. The top righthand panel shows the grayscale map of the 
DTFE density field reconstruction in and around the void, the corresponding velocity vector 
plot is shown in the bottom lefthand panel. Notice the detailed view of the velocity field: 
within the almost spherical global outflow of the void features can be recognized that can be 
identified with the diluted substructure within the void. Along the solid line in these panels we 
determined the linear DTFE density and velocity profile (bottom righthand frame). We can 
recognize the global "bucket" shaped density profile of the void, be it marked by substantial 
density enhancements. The velocity field reflects the density profile in detail, dominated by a 
global super-Hubble outflow. From Schaap 2007. 
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the void is dominated by the outflow of matter from the void, culminating into 
the void's own expansion near the outer edge. The comparison with the top two 
frames demonstrates the strong relation with features in the particle distribution and 
the DTFE density field. Not only is it slightly elongated along the direction of the 
void's shape, but it is also sensitive to some prominent internal features of the void. 
Towards the "SE" direction the flow appears to slow down near a ridge, while near 
the centre the DTFE reconstruction identifies two expansion centres. 

The general characteristics of the void expansion are most evident by following 
the density and velocity profile along a one-dimensional section through the void. 
The bottom-left frame of fig. [38] shows these profiles for the linear section along 
the solid line indicated in the other three frames. The first impression is that of 
the bucket-like shape of the void, be it interspersed by a rather pronounced density 
enhancement near its centre. This profile shape does confirm to the general trend 
of low-density regions to develop a near uniform interior density surrounded by 
sharply defined boundaries. Because initially emptier inner regions expand faster 
than the denser outer layers the matter distribution gets evened out while the inner 
layers catch up with the outer ones. The figure forms a telling confirmation of DTFE 
being able to recover this theoretically expected profile by interpolating its density 
estimates across the particle diluted void (see sect. 113. il l. 

The velocity profile strongly follows the density structure of the void. The lin- 
ear velocity increase is a manifestation of its general expansion. The near constant 
velocity divergence within the void conforms to the super-Hubble flow expected 
for the near uniform interior density distribution (see sect. ll3.lT ). Because voids are 
emptier than the rest of the universe they will expand faster than the rest of the 
universe with a net velocity divergence equal to 

= ^ = 3(a- 1), = i/void///, (48) 

ti 

where a is defined to be the ratio of the super-Hubble expansion rate of the void and 
the Hubble expansion of the universe. 

Evidently, the highest expansion ratio is that for voids which are completely 
empty, ie. /(void = -1- The expansion ratio a for such voids may be inferred from 
Birkhoff 's theorem, treating these voids as empty FRW universes whose expansion 
time is equal to the cosmic time. For a matter-dominated Universe with zero cos- 
mological constant, the maximum expansion rate that a void may achieve is given 
by 

= 1.5^3!;^ (49) 

with Qm the cosmological mass density parameter For empty voids in a Universe 
with a cosmological constant a similar expression holds, be it that the value of a 
will have to be numerically calculated from the corresponding equation. In gen- 
eral the dependence on A is only weak. Generic voids will not be entirely empty, 
their density deficit |zfvoidl ~ 0.8 - 0.9 (cf. eg. the linear density profile in fig. [38] ). 
The expansion rate 9„,id for such a void follows from numerical evaluation of the 
expression 
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(1 + 



(50) 



in which Qm,mid is the effective cosmic density parameter inside the void. 
13.2 Velocity Field Evolution 

On a mildly nonlinear scale the velocity field contains important information on 
the cosmic structure formation process. Density and velocity perturbations are sup- 
posed to grow gravitationally out of an initial field of Gaussian density and velocity 
perturbations. Once the fluctuations start to acquire values in the order of unity or 
higher the growth rapidly becomes more complex. The larger gravitational acceler- 
ation exerted by the more massive structures in the density field induce the infall of 
proportionally large amounts of matter and hence to an increase of the growth rate, 
while the opposite occurs in and around the underdense regions. 

The overdensities collapse into compact massive objects whose density excess 
may achieve values exceeding unity by many orders of magnitude. Meanwhile the 
underdensities expand and occupy a growing fraction of space while evolving into 
empty troughs with a density deficit tending towards a limiting value of -1.0. The 
migration flows which are induced by the evolving matter distribution are evidently 
expected to strongly reflect the structure formation process. 

In practice a sufliciently accurate analysis of the nonlinear cosmic velocities 
is far from trivial. The existing estimates of the peculiar velocities of galaxies are 
still ridden by major uncertainties and errors. This is aggravated by the fact that the 
sampling of the velocity field is discrete, rather poor and diluted and highly inhomo- 
geneous. While the conventional approach involves a smoothing of the velocity field 
over scales larger than lO/i 'Mpc in order to attain a map of the linear velocity field, 
we have argued in section [TO.41 that implicitly this usually yields a mass-weighted 
velocity field. DTFE represents a major improvement on this. Not only does it off'er 
an interpolation which guarantees an optimal resolution, thereby opening the win- 
dow onto the nonhnear regime, it also guarantees a volume-weighted flow field (see 
sect.[36]l. 

Density and Velocity Divergence 

The implied link between the matter distribution and the induced migration flows is 
most strongly expressed through the relation between the density field 6(x and the 
velocity divergence field. The bottom frames of fig. [39] contain greyscale maps of 
the DTFE normalized velocity divergence estimate 6 (with Hq the Hubble constant). 



e 



V ■ V _ 1 ( dvx dvy dv^ 



(51) 



Ho H[) [ dx dy 5z j ' 



for an N-body simulation. For illustrative purposes we have depicted the regions of 
negative and positive velocity divergence separately. The comparison between these 
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Fig. 39. The DTFE reconstructed velocity and velocity divergence field of a (low-resolution) 
SCDM N-body simulation, compared with the corresponding DTFE density field (top right- 
hand frame). The velocity divergence field is split in two parts. The negative velocity diver- 
gence regions, those marked by inflow, are shown in the bottom lefthand frame. They mark 
both the high-density cluster regions as well as the more moderate filamentary regions. The 
positive divergence outflow regions in the bottom righthand panel not only assume a larger 
region of space but also have a more roundish morphology. They center on the large (expand- 
ing) void regions in the matter distribution. The figure represents a nice illustration of DTFE's 
ability to succesfully render the non-Gaussian nonlinear velocity field. From Romano-Diaz 
2004. 



maps and the corresponding density field (upper righthand frame) and velocity field 
(upper lefthand) provides a direct impression of their intimate relationship. The two 
bottom frames clearly delineate the expanding and contracting modes of the veloc- 
ity field. The regions with a negative divergence are contracting, matter is falling in 
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along one or more directions. The inflow region delineates elongated regions coin- 
ciding with filaments and peaks in the density field. The highest inflow rates are seen 
around the most massive matter concentrations. Meanwhile the expanding regions 
may be seen to coincide with the surrounding large and deep underdense voids, 
clearly occupying a larger fraction of the cosmic volume. 

The strong relation between the density contrast and velocity divergence is a 
manifestation of the continuity equation. In the hnear regime this is a strict linear 
one-to-one relation. 



V ■ v(x, t) 
H 



= -a{t)f{Q„„Q^) 6{x,t), 



(52) 



linking the density perturb ation field 6 to the peculiar velocity field v via the factor 
f(Qm) (see iPeeble l ll980l) . There remams a 1-1 relation between velocity diver- 



gence and density into the mildly nonlinear regime (see eqn.l53]l. This explains why 
the map of the velocity divergence in fig.[39]is an almost near perfect negative image 
of the density map. 

Even in the quasi-linear and mildly nonlinear regime the one-to-one correspon- 
dance between velocity divergence and densjty remains intact, be it that it involves 
higher order terms (see iBernardeau et al.. 2002L for an exte nsive review). Within the 
context of Eulerian perturbation theory iBernardeaul d 1 992l (B)) derived an accurate 
2nd order approximati on form t he rel ation between the divergence and the density 
perturbation 5{x). iNusser et al.l (ll99U (N)) derived a similar quasi-linear approxi- 
mation within the context of the Lagrangian Zel'dovich approximation. According 
to these approximate nonlinear relations. 



1 

H 



V • v(x) = 



'lm„)[\-(i + 5(iL)fi^] (B) 

~f^^'"^ 1 +ai8 5(x) ^^"^ 



(53) 



for a Uni verse with Hubble parameter H(t) and r natter density parameter Qm. The 
studie s by^Bernardeau & van de Wevgaerl ( 19961) : Romano-Dfaz & van de Wevgaert 

;[sc 



(l2007h : lSchaap (2007.) have shown that the DTFE velocity field reconstructions are 
indeed able to reproduce these quasi-linear relations. 



Velocity Field Statistics 



The asymmetric nonlinear evolution of the cosmic velocity and density field mani- 
fests itself in an increasing shift of the statistical distribution of density and velocity 
perturbations away from the initial Gaussian probability distribution function. The 
analytical framework of Eulerian perturbation theor y provides a reasonably accu- 



rate description for the early nonlinear evolution (see iBernardeau et al.L 120021 for a 



review). As for the velocity field, robust results on the evolution and distribution of 
the velocity divergence, V ■ v, were derived in a series of papers by Bernardeau (e.g. 



Bernardeau et al.Lll995h . The complete probability distribution function (pdf) of the 
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velocity divergence may be evaluated via the summation of the series of cumulants 
of the non-Gaussian distribution function. The velocity divergence pdf is strongly 
sensitive to the underlying cosmological parameters, in particular the cosmological 
density parameter It represents a considerable advantage to the more regular 
analysis of the cosmic velocity flows on large linear scales. While the latter yields 
constraints on a combined function of Q„i and the bias b between the galaxy and 
matter distribution, access to nonlinear scales would break this degeneracy. 

An impression of the generic behaviour and changing global shape of the re- 
sulting non-Gaussian pdf as a function of Q,„ and the amplitude erg of the velocity 
divergence perturbations may be obtained from the top frame of fig.HO] The curves 
correspond to the theoretically predicted velocity divergence pdf for three differ- 
ent (matter-dominated) FRW cosmologies. Not only does Q„, influence the overall 
shape of the pdf, it also changes the location of the peak - indicated by 9ms.x - as 
well as that of the cutoff" at the high positive values of 9. By means of the Edgeworth 
expansion one m ay show that the velocity divergenc pdf reaches its ma ximum for a 



peak value Speak dJuszkiewicz et al.lll995l : lBernardeau & Kofmani[l995i) . 



epeak = -^cre; {O") ^ T^iO'f . (54) 
where the co efficient relates the 3rd order moment of the pdf to the second mo- 



ment (see e . g iBernardeaul 1 1 9941) . While the exact location of the peak is set by the 
bias-independent coefficient Tt„ it does not only include a dependence on Q„„ but 
also on the shape of the power spectrum, the geometry of the window function that 
has been used to filter the velocity field and on the value of the cosmological con- 
stant A. To infer the value of extra assumptions need to be invoked. Most direct 
therefore is the determination of Q,„ via the sharp cutoff of the pdf, related to the 
expansion velocity of the deepest voids in a particular cosmology (see eqn.|49]l. 

While conventional techniques may reproduce the pdf for moderate values of 
the velocity divergence 9, they tend to have severe difficulties in tracing the distri- 
bution for the more extreme positive values and the highest inffow velocities. An 
illustration of the tendency to deviate from the analytically predicted distribution 
can be seen in the two lefthand frames of fig. |40l showing the velocity divergence 
pdf determined from a SCDM N-body simulation for two scales of tophat filtered 
velocity fields {R = 10/z"'Mpc and R - 15/2 'Mpc). The open squares depict the 
velocity divergence pdf determined from an N-body simulation through a two-step 
grid procedure (see sect. 110. 4T i. Conventional grid interpolation clearly fails by huge 
margins to recover inffow onto the high-density filamentary structures as well as the 
outffow from voids. 

What came as the first evidence for the success of tessellation based interpola- 
tion is the rather tight correspondence of the Delaunay tessellation interpolation re- 



sults w ith the analytically predicted pdf. This finding bv lBernardeau & van de Wevgaerl 



19961) suggested that it would indeed be feasible to probe the non linear velocity 



field a nd infer directly accur4ate estimates of In a follow-up studv lBernardeau et al 



(Il997h succesfully tested this on a range of N-body simulations of structure forma- 



tion, showing Delaunay interpolation indeed recovered the right values for Q,„. The 
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Fig. 40. DTFE determination of the probability distribution function (pdf) of tfie velocity di- 
vergence 6. Top frame illustrates the sensitivity of the theta pdf to the underlying cosmology. 
The superposition of the theoretical pdf curves of three cosmologies immediately shows that 
the location of the maximum of the pdf and its sharp cutoff at the positive end are highly 
sensitive to Q. The other four frames demonstrate DTFE's ability to succesfuUy trace these 
marks. Lefthand frames: confrontation of the DTFE velocity divergence pdf and that deter- 
mined by more conventional two-step fixed grid interpolation method. Both curves concern 
the same 128^ SCDM N-body particle simulation (i.e. Q = 1). Top lefthand frame: tophat 
filter radius Rjn = lO/i 'Mpc. Bottom lefthand panel: R-^ = 15/7"'Mpc. The solid lines rep- 
resent theoretical predictions of the PDF for the measured values of erg (Bemardeau 1994). 
Black triangles are the pdf values measured by the DTFE method, the black squares by the 
equivalent VTFE Voronoi method. The open squares show the results obtained by a more 
conventional two-step fixed grid interpolation method. Especially for low and high 9 val- 
ues the tessellation methods prove to produce a much more genuine velocity divergence pdf. 
From Bemardeau & van de Weygaert 1996. Righthand frames: on the basis of DTFE's ability 
to trace the velocity divergence pdf we have tested its operation in a Q = 1 and a Q = 0.4 
CDM N-body simulation. For both configurations DTFE manages to trace the pdf both at its 
high- value edge and near its peak. From Bemardeau et al. 1997 



centre and bottom righthand frames depict two examples: the pdf 's for n Q — I and 
aQ - OA Universe accurately traced by the Delaunay interpolated velocity field. 



13.3 PSCz velocity field 



In a recent study iRomano-Dfaz & van de Wevsaertl ( 120070 applied DTFE towards 
reconstructing the velocity flow map throughout the nearby Universe v olume sam- 
pled by the PSCz galaxy redshift survey (also see iRomano-Diazl 120041) . 



The PSCz sample of the Local Universe 



The IRAS-PSCz catalo g JSaunders et al.Ll2000l) is an extension of the 1 .2-Jy cata- 
log ( IFisher et al.L Il995h . It contains ~ 15 500 galaxies with a flux at 60//m larger 



than 0.6-Jy. For a full description of the catalog, selection criteria and the proce- 
du res used to avo i d stell ar contamination and galactic cirrus, we refer the reader 
to Saunders et al. I (l2000l) . For our purposes the most important characteristics of 
the catalog are the large area sampled (~ 84% of the sky), its depth with a me- 
dian redshift of 8 500 km s"', and the dense sampling (the mean galaxy separation 
at 10 000 kms"' is </) - 1 000 kms"'). It impHes that PSCz contains most of 
the gravitationally relevant mass concentrations in our local cosmic neighbourhood, 
certainly sufficient to explain and study the cosmic flows within a local volume of 
radius ~ 120/7"'Mpc. 

To translate the redshift space distribution of galaxies in the PSCz catalog into 
galaxy positions and velocities i n real space the study b ased itself on an iterative pro- 
cessing of the galaxy sample by Branc hini et alj (1 19991) based on the linear theory of 
gravitational instability (Peebles^, 19801) . The method involved a specific realization 
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of an iterative technique to minimize redshift-space distortions dYahil et al.lll99l[) . 



While the resuhing galaxy velocities relate by construction to the linear clustering 
regime, the reconstruction of the velocity field throughout the sample volume does 
appeal to the capability of DTFE to process a discretely sampled field into a contin- 
uous volume-filling field and its ability to resolve flows in both high-density regions 
as well as the outflows from underdense voids. 



The PSCz velocity and density field 

The central circular field of fig. |4T] presents the DTFE velocity field in the Super- 
galactic Plane. For comparison, the bottom lefthand insert shows the discrete sam- 
ple of galaxy velocities which formed the input for the reconstruction. The velocity 
field is shown by means of the projected velocity vectors within the Z-supergalactic 
plane, superposed upon the corresponding DTFE density contourmaps inferred from 
the same PSCz galaxy sample. The length of the velocity arrows can be inferred 
from the arrow in the lower lefthand corner, which corresponds to a velocity of 650 
km/s. The Local Group is located at the origin of the map. 

The map of fig. |4T] shows the success of DTFE in converting a sample of dis- 
cretely sampled galaxy velocities and locations into a sensible volume-covering flow 
and density field. The processed DTFE velocity field reveals intricate details along 
the whole volume. The first impression is that of the meticulously detailed DTFE 
flow field, marked by sharply defined flow regions over the whole supergalactic 
plane. Large scale bulk flows, distorted flow patterns such as shear, expansion and 
contraction modes of the velocity field are clear features revealed through by the 
DTFE technique. DTFE recovers clearly outlined patches marked by strong bulk 
flows, regions with characteristic shear flow patterns around anisotropically shaped 
supercluster complexes, radial inflow towards a few massive clusters and, perhaps 
most outstanding, strong radial outflows from the underdense void regions. 

Overall, there is a tight correspondence with the large scale structures in the 
underlying density distribution. While the density field shows features down to a 
scale of Vs/i^'Mpc, the patterns in the flow field clearly have a substantially larger 
coherence scale, nearly all in excess of lOh'^Mpc. The DTFE velocity flow sharply 
follows the elongated ridge of the Pisces-Perseus supercluster. In addition we find 
the DTFE velocity field to contain markedly sharp transition regions between void 
expansion and the flows along the body of a supercluster. 

The local nature of the DTFE interpolation guarantees a highly resolved flow 
field in high density regions. Massive bulk motions are concentrated near and around 
the massive structure extending from the Local Supercluster (center map) towards 
the Great Attractor region and the Shapley concentration. The DTFE map nicely 
renders this pronounced bulk flow towards the Hydra-Centaurus region and shows 
that it dominates the general motions at our Local Group and Local Supercluster 
The top righthand insert zooms in on the flow in and around this region. The most 
massive and coherent bulk flows in the supergalactic plane appear to be connected 
to the Sculptor void and the connected void regions (towards the lefthand side of 



Geometric Analysis Cosmic Web 



93 



Fig. 41. Density and velocity field in the local Universe determined by DTFE on the basis 
of the PSCz galaxy redshift survey. Our Local Group is at the centre of the map. To the 
left we see the Great Attractor region extending out towards the Shapley supercluster. To 
the lefthand side we can find the Pisces-Perseus supercluster. The peculiar velocities of the 
galaxies in the PSCz galaxy redshift catalogue were determined by means of a linearization 
procedure (Branchini et al. 1999). The resulting galaxy positions and velocities (vectors) of 
the input sample for the DTFE reconstruction are shown in the bottom lefthand frame. The 

density values range from ~ 4.9 (red) down to 0.75 (darkblue), with cyan coloured regions 

having a density near the global cosmic average (6 ~ 0). The velocity vectors are scaled such 
that a vector with a length of x I /33rd of the region's diameter corresponds to 650 km/s. The 
density and velocity field have an effective Gaussian smoothing radius of Ra ~ VS/r'Mpc. 
The top righthand insert zooms in on the Local Supercluster and Great Attractor complex. 
From: Romano-Dfaz 2004. 



the figure). They are the manifestation of the combination of gravitational attrac- 
tion by the heavy matter concentration of the Pavo-Indus-Telescopium complex 
(at [SGX,SGY];^ [-40, -10]/! 'Mpc), the more distant "Hydra-Centaurus-Shapley 
ridge", and the effective push by the Sculptor void region. Conspicuous shear flows 
can be recognized along the ridge defined by the Cetus wall towards the Pisces- 
Perseus supercluster ([SGX,SGY]a; [20, -40]/!" 'Mpc. A similar strong shear flow 
is seen along the extension of the Hydra-Centaurus supercluster towards the Shapley 
concentration. 

The influence of the Coma cluster is beautifully outlined by the strong and near 
perfect radial infall of the surrounding matter, visible at the top-centre of figure [41] 
Also the velocity field near the Perseus cluster, in the Pisces-Perseus supercluster 
region, does contain a strong radial inflow component. 

Perhaps most outstanding are the radial outflow patterns in and around voids. 
In particular its ability to interpolate over the low-density and thus sparsely sam- 
pled regions is striking: the voids show up as regions marked by a near-spherical 
outflow. The intrinsic suppression of shot noise effects through the adaptive spa- 
tial interpolation procedure of DTFE highlights these important components of the 
Megaparsec flow field and emphasizes the dynamical role of voids in organizing the 
matter distribution of the large scal e Cosmic Web. By contrast, more c onventional 
schemes, such as TSC or SPH (see ISchaap & van de Wevsaertl l2007al) . meet sub- 
stantial problems in defining a sensible field reconstruction in low density regions 
without excessive smoothing and thus loss of resolution. 



14 DTFE meets reality: 2dFGRS and the Cosmic Web 

To finish the expose on the potential of the Delaunay Tessellation Field Estimator, 
we present a reconstruction of the foamy morphology of the galaxy distribution 
in the 2dF Galaxy Redshift Survey (2dFGRS). DTFE was used to reconstruct the 
projected galaxy surface density field as well as the full three-dimensional galaxy 
density field. 
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Fig. 42. DTFE galaxy surface density reconstructions of tlie projected galaxy distribution in 
the two 2dF slices (north and south). For comparison see the galaxy distribution in fig. [T] A 
description may be found in the text (sect. I14.2t . From Schaap 2007 and Schaap & van de 
Weygaert 2007d. 



Geometric Analysis Cosmic Web 



95 



14.1 the 2dF Galaxy Redshift Survey 

The 2dFGRS is one of the major spectroscopic surveys in which the spectra of 
245 591 objects have been obtained, with the scope of obtaining a repr esentative pic 



tureof the large scale distribution of galaxies in the nearby universe (IColless et al 



20031) . It is a magnitude-limited survey, with galaxies selected do wn to a limiting 



magni tude of bj ~ 19.45 from the extended APM Galaxy Survey (IMaddox et al 
1990allblH) . The galaxy redshifts of galaxies were measured with the 2dF multifibre 
spectrograph on the Anglo- Australian Telescope, capable of simultaneously observ- 
ing 400 objects in a 2° diameter field. The galaxies were confined to three regions, 
together covering an area of approximately 1500 squared degrees. These regions 
include two declination strips, each consisting of overlapping 2° fields, as well as 
a number of 'randomly' distributed 2° control fields. One strip (the SOP strip) is 
located in the southern Galactic hemisphere and covers about 80° x 15° close to the 
South Galactic Pole (2l'W" <a< 03'W", -37.5° <6< -22.5°). The other strip 
(the NGP strip) is located in the northern Galactic hemisphere and covers about 
75° X 10°(09''50'" <a< 14*50"', -7.5° <6 < +2.5°). Reliable redshifts were ob- 
tained for 221 414 galaxies. 

These data have been made publically available in the form of the 2dFGRS final 
data release (available at |http : //msowww . anu . edu . au/2dFGRS/). 

14.2 Galaxy surface density field reconstructions 

The galaxy distribution in the 2dFGRS is mainly confined to the two large strips, 
NGP and SGP Since their width is reasonably thin, a good impression of the spatial 
patterns in the galaxy distribution may be obtained from the 2-D projection shown 
in fig. m 

We have reconstructed the galaxy surface density field in redshift space in the 
2dFGRS NGP and SGP regions. All density field reconstructions are DTFE recon- 
structions on the basis of the measured galaxy redshift space positions. As no cor- 
rections were attempted to translate these into genuine positions in physical space, 
the density reconstructions in this section concern redshift space. In order to war- 
rant a direct comparison with the galaxy distribution in fig. [T] the results shown in 
this section were not corrected for any observational selection eff'ect, also not for the 
survey radial selection fun ction. For selec t ion corrected density field reconstru ctions 



we refer to the analysis in lSchaapl (l2007h : ISchaap & van de Weygaertj (l2007cl) 



Fig. |42] shows the resulting DTFE reconstructed density field. DTFE manages 
to reveal the strong density contrasts in the large scale density distribution. The 
resolution is optimal in that the smallest interpolation units are also the smallest 
units set by the data. At the same time the DTFE manages to bring out the fine 
structural detail of the intricate and often tenuous filamentary structures. Particularly 
noteworthy are the thin sharp edges surrounding voidlike regions. 
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Fig. 43. DTFE galaxy surface density in selected regions in the 2dFGRS galaxy surface den- 
sity field. For the density field in the total 2dFGRS region see fig. |42] For a discussion see 
sect. 114.21 Frame 1 zooms in on the Great Wall in the southern (SGP) 2dF slice, frame 5 on 
the one in the northern (NGP) slice. Note the sharply rendered "fingers of God" marking the 
sites of clusters of galaxies. Perhaps the most salient feature is the one seen in frame 3, the 
cross-like configuration in the lower half of the NGP region. From Schaap 2007 and Schaap 
& van de Weygaert 2007d. 



Individual Structures 

The impressive resolution and shape sensitivity of the DTFE reconstruction be- 
comes particularly visible when zooming in on structural details in the Cosmic Web. 
Figure |43] zooms in on a few interesting regions in the map of fig.[T] Region 1 fo- 
cuses on the major mass concentration in the NGP region, the Sloan Great Wall. 
Various filamentary regions emanate from the high density core. In region 2 a void- 
like region is depicted. The DTFE renders the void as a low density area surrounded 
by various filamentary and wall-like features. Two fingers of God are visible in the 
upper right-hand corner of region 2, which show up as very sharply defined lin- 
ear features. Many such features can be recognized in high density environments. 
Note that the void is not a totally empty or featureless region. The void is marked 
by substructure and contai ns a number of smaller subvoids, reminding of the hier - 



archical evolution of voids jDubinski et a/.LI1993l:lSheth & van de Wevgaern, 120041) 



Region 3 is amongst the most conspicuous structures encountered in the 2dF survey. 
The cross-shaped feature consists of four tenuous filamentary structures emanating 
from a high density core located at the center of the region. Region 4 zooms in on 
some of the larger structures in the SGP region. Part of the Pisces-Cetus superclus- 
ter is visible near the bottom of the frame, while the concentration at the top of this 
region is the upper part of the Horologium-Reticulum supercluster. Finally, region 5 
zooms in on the largest structure in the SGP region, the Sculptor supercluster. 



DTFE artefacts 

Even though the DTFE offers a sharp image of the cosmic web, the reconstructions 
also highlight some artefacts. At the highest resolution we can directly discern the 
triangular imprint of the DTFE kernel. Also a considerable amount of noise is visible 
in the reconstructions. This is a direct consequence of the high resolution of the 
DTFE reconstruction. Part of this noise is natural, a result of the statistical nature of 
the galaxy formation process. An additional source of noise is due to the fact that the 
galaxy positions have been projected onto a two-dimensional plane. Because DTFE 
connects galaxies which lie closely together in the projection, it may involve objects 
which in reality are quite distant. A full three-dimensional reconstruction followed 
by projection or smoothing on a sufficiently large scale would alleviate the problem. 
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Fig. 44. Isodensity surface of the galaxy distribution in the north (top) and south region 
(bottom) of the 2dFGRs. The densities are determined by means of the DTFE technology, 
subsequently Gaussian smoothed on a scale of 2/r'Mpc. Several well-known structures are 
indicated. From Schaap 2007 and Schaap & van de Weygaert 2007d. 
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14.3 Three-dimensional structure 2dFGRS 

We have also reconstructed the full three-dimensional galaxy density field in the 
NGP and SGP regions of the 2dFGRS. The result is shown in Fig. IH It shows 
the three-dimensional rendering of the NGP (left) and SGP slices (right) out to a 
redshift z - 0.1. The maps demonstrate that DTFE is able of recovering the three- 
dimensional structure of the cosmic web as well as its individual elements. Although 
less obvious than for the two-dimensional reconstructions, the effective resolution 
of the three-dimensional reconstructions is also varying across the map. 

The NGP region is dominated by the large supercluster at a redshift of about 
0.08, the Sloan Great Wall. The structure near the upper edge at a redshift of 0.05 to 
0.06 is part of the upper edge of the Shapley supercluster In the SGP region several 
known superclusters can be recognized as well. The supercluser in the center of this 
region is part of the Pisces-Cetus supercluster The huge concentration at a redshift 
of about 0.07 is the upper part of the enormous Horologium-Reticulum supercluster 

15 Extensions, Applications and Prospects 

In this study we have demonstrated that DTFE density and velocity fields are op- 
timal in the sense of defining a continuous and unbiased representation of the data 
while retaining all information available in the point sample. 

In the present review we have emphasized the prospects for the analysis of we- 
blike structures or, more general, any discretely sampled complex spatial pattern. 
In the meantime, DTFE has been applied in a number of studies of cosmic struc- 
ture formation. These studies do indeed suggest a major improvement over the more 
conventional analysis tools. Evidently, even though density/intensity field analysis 
is one of the primary objectives of the DTFE formalism one of its important fea- 
tures is its ability to extend its Delaunay tessellation based spatial interpolation to 
any corresponding spatially sampled physical quantity. The true potential of DTFE 
and related adaptive random tessellation based techniques will become more ap- 
parent as further applications and extensions will come to the fore. The prospects 
for major developments based on the use of tessellations are tremendous. As yet 
we may identify a diversity of astrophysical and cosmological applications which 
will benefit substantially from the availability of adaptive random tessellation re- 
lated methods. A variety of recent techniques have recognized the high dynamic 
range and adaptivity of tessellations to the spatial and morphological resolution of 
the systems they seek to analyze. 

The first major application of DTFE concerns its potential towards uncover- 
ing morphological and statistical characteristics of spatial patterns. A second major 
class of applications is that of issues concerning the dynamics of many particle sys- 
tems. They may prove to be highly useful for the analysis of phase space structure 
and evolution of gravitational systems in their ability to efliciently trace density 
concentrations and singularities in higher dimensional space. Along similar lines a 
highly promising avenue is that of the application of similar formalisms within the 
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context of computational hydrodynamics. The apphcation of Delaunay tessellations 
as a fully self-adaptive grid may finally open a Lagrangian grid treatment of the 
hydrodynamical equations describing complex multiscale systems such as encoun- 
tered in cosmology and galaxy formation. Various attempts along these lines have 
already been followed and with the arrival of efficient adaptive Delaunay tessella- 
tion calculations they may finally enable a practically feasible implementation. A 
third and highly innovative application, also using both the adaptive and random 
nature of Delaunay tessellations, is their use as Monte Carlo solvers for systems 
of complex equations descri bing complex phys ical systems. Particularly interesting 
has been the recent work by Ritzerveld & Ickel ((2006). They managed to exploit the 
random and adaptive nature of Delaunay tessellations as a stochastic grid for Monte 
Carlo lattice treatment of radiative transfer in the case of multiple sources. In the 
following we will describe a few of these applications in some detail. 



15.1 Gravitational Dynamics 



Extrapolating the observatio n that DTFE is able to simultaneously handle spatial 
density a nd velocity field (e.g .^Bernardeau & va n de Weveaerl , 1996 ; Romano-Dfaz , 

12004: ,Romano-Dfaz & van de Weygaert. ,2007.) . and encouraged by th e success of 

Voronoi-based methods i n identifying dark halos in N-body simulations ( Neyrinck. Gnedin & Hamiltoni 



20051) . lArad et al.l (120041) used DTFE to assess the six-dimensional phase space den- 



sity distribution of dark halos in cosmological A^-body simul ations. While a fully 
six-di mensional analysis may be computationally cumbersome dAscasibar & Binnev , 



2005|), and not warranted because of the symplectic character of phase-space, the 



splitting of the problem into a separate spatial and velocity-space three-dimensional 
tessellation may indeed hold promise for an innovative analysis of the dynamical 
evolution of dark halos. 



Gravitational Lensing and Singularity Detection 



A related promising ave nue seeks to use the ability of DTFE to trace sharp density 
contrasts. This impelled iBradac et al.l ( 120041) to apply DTFE to gravitational lens- 
ing. They computed the surface density map for a galaxy from the projection of the 
DTFE volume density field and used the obtained surface de nsity map to co mpute 
the gravitational lensing pattern around the object. Recently, iLi et al. [ (l2006l) eval- 
uated the method and demonstrated that it is indeed a promising tool for tracing 
higher-order singularities. 



15.2 Computational Hydrodynamics 

Ultimately, the ideal fluid dynamics code would combine the advantages of the Eu- 
lerian as well as of the Lagrangian approach. In their simplest formulation, Eulerian 
algorithms cover the volume of study with a fixed grid and compute the fluid trans- 
fer through the faces of the (fixed) grid cell volumes to follow the evolution of the 
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system. Lagrangian formulations, on the other hand, compute the system by follow- 
ing the ever changing volume and shape of a particular individual element of gafl 
Their emphasis on mass resolution makes Lagrangian codes usually better equipped 
to follow the system into the highest density regions, at the price of a decreased res- 
olution in regions of a lower density. 

As we may have appreciated from the adaptive resolution and morphological 
properties of Delaunay tessellations in DTFE and the more general class of Natural 
Neighbour interpolations the recruitment of Delaunay tessellations may define an 
appropriate combination of Eulerian and Lagrangian description of fluid dynamical 
systems. 



Particle Hydrodynamics 



A well-known strategy in computational hydrodynamics is to follow the Lagrangian 
description by discretizing the fluid system in terms of a (large) number of fluid par- 
ticles which encapsulate, trace and sample t he relevant dynamical and thermody- 
namical characteristics of the entire fluid (see Koumoutsakosl 2005 , for a recent re- 
view) . Smooth Particle H ydrodynamics (SPH) codes ( ILucv , 19771 : Gingold & Monaghan , 
1977 ; Monaghan! 1992 ) have found widespread use in many areas of science and 
have arguably become the most prominent computational tool in cosmology (e.g. 



Katz. Weinberg & HernquisIll996l:[Springel. Yoshida & WhiteLEoOltlSpringelLbOOsI) . 
SPH particles should be seen as discrete balls of fluid, whose shape, extent and prop- 
erties is specified according to user-defined criteria deemed appropriate for the phys- 
ical system at hand. Notwithstanding substantial virtues of SPH - amongst which 
one should include its ability to deal with systems of high dynamic range, its adap- 
tive spatial resolution, and its flexibility and conceptual simplicity - it also involves 
various disadvantages. 

One straightforward downside concerns its comparatively bad performance in 
low-density regions. In a highly topical cosmological issue such as that of the reion- 
ization of the Universe upon the formation of the first galaxies and stars it may 
therefore not be able to appropriately resolve the void patches which may be rela- 
tively important for the transport of radiation. Even in high density regions it may 
encounter problems in tracing singularities. As a result of its user-defined finite size 
it may not succeed in a sufficiently accurate outlining of shocks. The role of the 
user-dependent and artificial particle kernel has also become apparent when com- 
paring the performance of SPH versus the more natural DTFE kernels (see fig. 112. II 
fig.IUJJand sect.innj. 

In summary: SPH involves at least four user-defined aspects which afl'ect its 
performance. 



1 . SPH needs an arbitrary user-specified kernel function W. 



^Interestingly, the Lagrangian formulation is also due to Euler l lEuleJ.ll862h who em- 
ployed this formalism in a le tter to Lagrange, who later proposed these ideas in a publication 
by himself l lLagrang3 , ll762l) 
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2. SPH needs a smoothing length h; 

even though the standard choice is a length adapting to the particle density it 
does imply a finite extent. 

3. SPH kernel needs a shape; a spherical shape is usually presumed. 

4. SPH needs an (unphysical) artificial viscosity to stabilize solutions and to be 
able to capture shocks. 

Given the evidently superior density and morphology tracing abilities of DTFE 



and re lated methods based upon adaptive grids lPelupessy. Schaap & van de Weygaerll 



(120031) investigated the question what the effect would be of replacing a regular SPH 
kernel by an equivalent DTFE based kernel. In an application to a TreeSPH simula- 
tion of the (neutral) ISM they managed to identify various density and morphologi- 
cal environments where the natural adaptivity of DTFE proves to yield substantially 
better results. They concluded with various suggestions for the formulation of an 
alternative particle based adaptive hydrodynamics code in which the kernel would 
be defined by DTFE. 

A closely related and considerably advanced development concerns the formu- 
lation of a complete and s elf-consistent particle hydrodynamics code by Espanol , 



panol, 

1 120001: 

Serrano & EspanoiboOltlSerrano et al.,.2002:.De Fabritiis & Coveney».2003D . Their 



Coveney and collaborators (lEspafio l. 1998; Flekk0v, Covenev & De Fabritiis 



(Multiscale) Dissipative Particle Hydrodynamics code is entirely formulated on the 
basis of Voronoi fluid particles. Working out the Lagrangian equations for a fluid 
system they demonstrate that the subsequent compartimentalization of the system 
into discrete thermodynamic systems - fluid particles - leads to a set of equations 
whose solution would benefit best if they are taken to be defined by the Voronoi 
cells of the corresponding Voronoi tessellation. In other words, the geometrical fea- 
tures of the Voronoi model is directly connected to the physics of the system by 
interpretation of the Voronoi volumes as coarse-grained "fluid clusters". Not only 
does their formalism capture the basic physical symmetries and conservation laws 
and reproduce the continuum for (sampled) smooth fields, but it also suggests a 
striking similarity with turbulence. Their formalism has been invoked in the mod- 
eling mesoscale systems, simulating the molecular dynamics of complex fluid and 
soft condense matter systems which are usually are marked fluctuations and Brow- 
nian motion. While the absence of long-range forces such as gravity simplifies the 
description to some extent, the Voronoi particle descriptions does provide enough 
incentive for looking into possible implementations within an astrophysical context. 

Adaptive Grid Hydrodynamics 

For a substantial part the success of the DTFE may be ascribed to the use of Delau- 
nay tessellations as optimally covering grid. This implies them to also be ideal for 
the use in moving &■ adaptive grid implementations of computational hydrodynam- 
ics. 

Unstructured grid methods originally emerged as a viable alternative to struc- 
tured or block-structured grid techniques for discretization of the Euler and Nav ier- 



Stokes equations for structures with complex geometries (see iMavripilisL 1 1 9971 for 
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Fig. 45. Application of the Natural Neighbour scheme solution of partial differential equa- 
tions on highly irregular evolving Delaunay grids, described by Braun & Sambridge 1995. 
It involves the natural-element method (NEM) solution of the Stokes flow problem in which 
motion is driven by the fall of an elasto-plastic plate denser than the viscous fluid. The prob- 
lem is solved by means of a natural neighbour scheme, the Delaunay grid is used as unstruc- 
tured computation mesh. The Stokes equation is solved at the integration points in the linear 
fluid, the equations of force balance at the integration points located within the plate. The 
solution is shown at two timesteps (1 and 10000 of the integration). Image courtesy of M. 
Sambridge also see Braun &. Sambridge 1995. 



an extensive review). The use of unstructured grids provides greater flexibility for 
discretizing systems with a complex shape and allows the adaptive update of mesh 
points and connectivities in order to increase the resolution. A notable and early 
example of the use of unstructu red grids may be found General Relativity in the 
formulation of Regge calculus (IReggel Il96lb . Its applications incl udes q i iantum 
gravity formulations on the basis of Voronoi and Delaunay grids (e.g. lMilleiill997b . 



One class of unstructured grids is based upon the use of specific geometrical 
elements (triangles in 2-D, tetrahedra in 3-D) to cover in a non-regular fashion a 
structure whose evolution one seeks to compute. It has become a prominent tech- 
nique in a large variety of technological applications, with those used in the design 
of aerodynamically optimally shaped cars and the design of aeroplanes are perhaps 
the most visible ones. The definition and design of optimally covering and suitable 
meshes is a major industry in computational science , including issues involved with 



the definition of optim al Delaunay meshe s (see e.glAmenta. B ern & K amvvsselis , 



19981; lAmeiita & Bernl 1 19991: IShewchtifl 120021; ICheng et al.i. .2004; .Devi 12004 . 
Alliez et al. , 2005allbt Boissonnat et al. , 20061) . A second class of unstructured grids 



is based upon the use of structural elements of a mixed type with irregular connec- 
tivity. The latter class allows the definition of a self-adaptive grid which follows the 
evolution of a physical system. It is in this context that one may propose the use of 
Voronoi and Delaunay tessellations as adaptive hydrodynamics lattices. 
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Hydrocodes with Delaunay tessellations at their core warrant a close connec- 
tion to the underlying matter distribution and would represent an ideal compromise 
between an Eulerian and Lagrangian description, guaranteeing an optimal resolu- 
tion and dynamic range while taking care of an appropriate coverage of low-density 
regions. Within astrophysical context there have been a few notable attempts to de 



velop moving grid codes. Although these have shown their potential (lGnedinill995 
|Pen, 1998), their complexity and the implicit complications raised by the dominant 
presence of the long-range force of gravity have as yet prevented their wide-range 
use. 

It is here that Delaunay tessellations may prove to offer a highly promising 
alternative. The advantages of a moving grid fluid dynamics code based on De- 
launay tessellations has been explicitly addressed in a detailed and enticing study 
by [whitehurst (1995). His two-dimensional FLAME Lagrangian hydrocode used 
a first-order solver and proven to be far superior to all tested first-order and many 
second-order Eulerian codes. The adaptive nature of the Langrangian method and 
the absence of preferr ed directions in the g rid proved to be key factors for the perfor- 
mance of the method. WhitehursS ( 1995 ) illustrated this with impressive examples 



such as the capturing of shock and collision of perfectl y spherical shells. T he related 
higher-order Natural Neighbour hydrocod es used by Braun & Sambrid ge ( 1995|, 



for a range of geophysical problems, and iSukumari (Il998h are perhaps the most 
convincing examples and applications of Delaunay grid based hydrocodes. The ad- 
vantages of Delaunay grids in principle apply to any algorithm invoking them, in 
particular also for three-dimensional implementations (of which we are currently 
unaware). 

Kinetic and Dynamic Delaunay Grids 

If anything, the major impediment towards the use of Delaunay grids in evolving 
(hydro)dynamical systems is the high computational expense for computing a De- 
launay grid. In the conventional situation one needs to completely upgrade a De- 
launay lattice at each timestep. It renders any Delaunay-based code prohibitively 
expensive. 

In recent years considerable attention has been devoted towards developing ki- 
netic and dynamic Delaunay codes which manage to limit the update of a grid 
to the few sample points that have moved so far that the identity of their Nat- 
ural Neighbours has changed. The work by Meyer-Hermann and collaborators 



dSchaller & Mever-Hermann 2004 : Bever et al. , 2005 : Bever & Mever-Hermann 
2006h may prove to represent a watershed in having managed to define a parallel 



code for the kinetic and dynamic calculation of Delaunay grids. Once the Delau- 
nay grid of the initial point configuration has been computed, subsequent timesteps 
involve a continuous dynamical upgrade via local Delaunay simplex upgrades as 
points move around and switch the identity of their natural neighbours. 

The code is an elaboration on the more basic step of inserting one point into 
a Delaunay grid and evaluating which Delaunay simplices need to be upgraded. 
Construction of Delaunay triangulations via incremental insertion was introduced 
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Fig. 46. Nine non-linear gray-scale images of the density evolution of the 2D interacting blast 
waves test of the FLAME Delaunay grid hydrocode of Whitehurst (1995). The circular shape 
of the shockfronts is well represented; the contact fronts are clearly visible. Of interest is the 
symmetry of the results, which is not enforced and so is genuine, and the instability of the 
decelerated contact fronts. From Whitehurst 1995. 



by iFortund (1 1 992h and lEdelsbrunner & Shahl(ll996h . Only sample points which ex- 
perience a change in identity of their natural neighbours need an update of their 
Delaunay simpUces, based upon vertex flips. The update of sample points that have 
retained the same natural neighbours is restricted to a mere recalculation of the po- 
sition and shape of their Delaunay cells. 
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15.3 Random Lattice Solvers and Simplex 

The use of Delaunay tessellations as adaptive and dynamic meshes for evaluating 
the fluid dynamical equations of an evolving physical system emphasizes their adap- 
tive nature. Perhaps even more enticing is the use of the random nature of the these 
meshes. Taking into account their spatial adaptivity as well as their intrinsic stochas- 
tic nature, they provide an excellent template for following physical systems along 
the lines of Monte Carlo methods. 




Fig. 47. The result of two simple Simplex radiative transfer tests on a 2D Poisson-Delaunay 
random lattice with N = 5x 10* points. Both are logarithmic plots of the number of particles 
at each site. Left: illustration of the conservation of momentum by means of the transport 
of particles with constant momentum vectors. Right: Illustration of a scattering transport 
process. Image courtesy of J. Ritzerveld, also see Ritzerveld 2007. 



Generally speaking, Monte Carlo methods determine the characteristics of many 
body systems as statistical averages over a large number of particle histories, which 
are computed with the use of random sampling techniques. They are particularly 
useful for transport problems since the transport of particles through a medium can 
be described stochastically as a random walk from one interaction event to the next. 
The first t o develop a computer-orien ted Monte Carlo method for transport prob- 
lems were [Metropolis & Ulam (Il949l) . Such transport problems may in general be 
formulated in terms of a stochastic Master equation which may evaluated by means 
of Monte Carlo met hods by simulatin g large numbers of different particle trajecto- 
ries or histories (see iRitzerveldl |2007[ for a nice discussion). 

While the asymmetry of regular meshes and lattices for Monte Carlo calcula- 
tions introduces various undesirable eff'ects, random lattices may alleviate or solve 
the problematic lack of symmet ry. Lattice Boltzmann s tudies were amongst the 
first to recognize this (see e.g. lUbertini & Succil l2005b . In a set of innovative 
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Fig. 48. A volume rendering of the result of using the Simplex method to transport ioniz- 
ing photons through the intergalactic medium at and around the epoch of reionization. The 
Simplex method was applied to a PMFAST simulation of the large scale structure in a ACDM 
universe. The results of 6 different timesteps are plotted, in which the white corresponds to 
hydrogen to the hydrogen that is still neutral (opaque), the blue to the already ionized hydro- 
gen (transparent). Image courtesy of J. Ritzerveld, also see Ritzerveld 2007. 
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publications Christ. Friedberg & Lee ( 1982albllc ) applied random lattices, including 
Voronoi and Delaunay tessellations, to solve (QCD) field theoretical equations. 

Recently, a highly innovative and interesting appUcation to (astrophysical) ra- 
diative transfer problems followed the sam e philosophy as t he Random Lattice 
Gauge theories. iRitzerveld & Icke (l2006h and lRitzerveldl ( 120071) translated the prob- 
lem of radiative transfer into one which would be solved on an appropriately de- 
fined (Poisson)-Delaunay grid. It leads to the formulation of the Simplex radiation 
transfer technique which translates the transport of radiation through a medium into 
the deterministic displacement of photons from one vertex to another vertex of the 
stochastic Delaunay mesh. The perfect random and isotropic nature of the Delaunay 
grid assures an unbiased Monte Carlo sampling of the photon paths. The specifica- 
tion of appropriate absorption and scattering coefficients at the nodes completes the 
proper Monte Carlo solution of the Master equation for radiative transfer 

One of the crucial aspects of Simplex is the sampling of the underlying density 
field p(x) such that the Delaunay grid edge lengths Lq represent a proper random 
sample of the (free) paths Ay of photons through a scattering and absorbing medium. 
In other words. Simplex is built on the assumption that the (local) mean Delaunay 
edge length (Ld) should be proportional to the (local) mean free path of the pho- 
ton. For a medium of density p(x) in d-dimensional space (d=2 or d=3), with ab- 
sorption/scattering coefficient cr, sampled by a Delaunay grid generated by a point 
sample with local number density «x(x). 



<Ld) = ^(^/)«x^' 



- p(x)<r 



(55) 



where ^(d) is a dimension dependent constant. By sampling of the underlying den- 
sity field p(x) by a point density 



(56) 



Simplex produces a Delaunay grid whose edges are guaranteed to form a represen- 
tative stochastic sample of the photon paths through the medium. 

To illustrate the operation of Simplex Fig. |47] shows the outcome of a two- 
dimensional Simplex test calculations. One involves a test of the ability of Simplex 
to follow a beam of radiation, the second one its ability to follow the spherical spread 
of radiation emitted by an isotropically emitting source. The figure nicely illustrates 
the success of Simplex, meanwhile providing an impression of the ebffect and its 
erratic action at the scale of Delaunay grid cells. 

The Simplex method has been applied to the challenging problem of the reion- 
ization of the intergalactic medium by the first galaxies and stars. A major prob- 
lem for most radiative transfer techniques is to deal with multiple sources of ra- 
diation. While astrophysical problems may often be a pproximated by a descrip- 
tion in terms of a single source, a proper evaluation of reionization should take 
into account the role of a multitude of sources. Simplex proved its ability to 
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yield sensible answers in a series of test comparisons bet ween different rad iative 
transfer codes applied to aspects typical for reionization (IllUev et al. , l2006h . The 
first results of the application of Simplex to genuine cosmological circumstances, 
by couphng it to a cosmolog i cal SP H code, do yield highly encouraging results 
dRitzerveld. Pawlik & Schave , l2007h . Fig. HHlis a beautiful illustration of the out- 
come of one of the reionization simulations. 

If anything. Simplex proves that the use of random Delaunay grids have the 
potential of representing a genuine breakthrough for addressing a number of highly 
complex astrophysical problems. 



15.4 Spatial Statistics and Pattern Analysis 



Within its cosmological context DTFE will meet its real potential in more sophis- 
ticated applications tuned towards uncovering morphological characteristics of the 
reconstructed spatial patterns. 

A potentially interesting application would be its implementation in the SURF- 
GEN machinery. SURFGEN seeks to provide locall y defined topological measures, 
cq. local Min kowski functionals, of the density fie ld (Sahn i. Sathvaprakash & Shandarin , 



ri998; Sheth et"all 120031: IShandarin. Sheth & Sahnl .20o4 . A recent new sophisti 



cat ed technique for tracing the cosm ic web is the skeleton formalism developed 
bv 'Novi kov. Colom bi & Pore ('2006'), based upon the framework of Morse theory 
fColombi. Pogosvan & Souradeep. 2OOO0 . The skeleton formalism seeks to identify 
filaments in the web by identifying ridges alo ng which the gradien t of the density 
field is extremal along its isocontours (also see lSousbie et al. [ l2006h . The use of un- 
biased weblike density fields traced by DTFE would undoubtedly sharpen the ability 
of the skeleton formalism to trace intricate weblike features over a wider dynamical 
range. 

Pursuing an alternative yet related track concerns the direct use of tessellations 
themselves in outlining topological properties of a complex spatial distribution of 
points. Alpha-shapes of a discrete point distribution are subsets of Delaunay triangu- 
lations which are sensitive tracers of its topology and may be exploited towards in- 
ferring Minkowski functionals and Betti numbers dEdelsbrunner. Kirkpatrick & Seidell 



1983tlEdelsbrunner & Muckelll994tlEdelsbrunner. Letscher & Zomorodianil2002h 
A rec ent and ongoing stu dy seeks their application to the description of the Cosmic 
Web dVegter et al.Ll2007h (also see sec ll5.5l i. 

Two major extensions of DTFE already set out to the identification of key as- 
pects of the cosmic web. Specifically focussed on the hierarchical nature of the 
Megaparsec matter distribution is the detection of weblike anisotropic features over 
a range of spa tial scales by the M ultiscale Morphology Filter (MMF), introduced 
and defined by Aragon-Calvo et al. ( 2007). The Wat ershed Void Finding (WVF) al- 
gorithm of IPlaten. van de Weygaert & JonesI (120071) is a void detection algorithm 
meant to outline the hierarchical nature of the cosmic void population. We will 
shortly touch upon these extensions in order to illustrate the potential of DTFE. 
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15.5 Alpha-shapes 



Alpha shape are a description of the (intuitive) notion of the shape of a dis- 
crete point set. Alpha-shapes of a discrete point distribution are subsets of a 
Delaunay triangulation and were introduced by Edelsbrunner and collaborators 
tede sbrunner. Kirkpatrick & Seidell 1 1 9831: iMuckel 1 1 9931: fedelsbrunner & Mucke . 



1994 



Edelsbrunner. Letscher & Zomorodiani l2002h . Alpha shapes are generaliza- 



tions of the convex hull of a point set and are concrete geometric objects which are 
uniquely defined for a particular point set. Reflecting the topological structure of a 
point distributio n, it is one of the most essential concepts in t he field of Computa- 



tional Topology jDeyi 119981: IVegteit 120041: IZomorodianl 120051) 



If we have a point set S and its corresponding Delaunay triangulation, we may 
identify all Delaunay simpUces - tetrahedra, triangles, edges, vertices - of the tri- 
angulation. For a given non-negative value of a, the alpha complex of a point sets 
consists of all simplices in the Delaunay triangulation which have an empty circum- 
sphere with squared radius less than or equal to a, < a. Here "empty" means that 
the open sphere does not include any points of S . For an extreme value a - Q the 
alpha complex merely consists of the vertices of the point set. The set also defines 
a maximum value am^x, such that for a > cfmax the alpha shape is the convex hull of 
the point set. 

The alpha shape is the union of all simplices of the alpha complex. Note that it 
implies that although the alpha shape is defined for all < ff < oo there are only 
a finited number of different alpha shapes for any one point set. The alpha shape is 
a polytope in a fairly general sense, it can be concave and even disconnected. Its 
components can be three-dimensional patches of tetrahedra, two-dimensional ones 
of triangles, one-dimensional strings of edges and even single points. The set of 
all real numbers a leads to a family of shapes capturing the intuitive notion of the 
overall versus fine shape of a point set. Starting from the convex hull of a point set 
and gradually decreasing a the shape of the point set gradually shrinks and starts to 
develop cavities. These cavities may join to form tunnels and voids. For sufficiently 
small a the alpha shape is empty. 

Following this description, one may find that alpha shapes are intimately related 
to the topology of a point set. As a result they form a direct and unique way of char- 
acterizing the topology of a point distribution. A complete quantitative description 
of the topology is that in terms of Betti numbers ySk and these may indeed be di- 
rectly inferred from the alpha shape. The first Betti number /3o specifies the number 
of independent components of an object. In R /3i may be interpreted as the number 
of independent tunnels, and ^62 as the number of independent enclose voids. The k''' 
Betti number effectively counts the number of independent ^-dimensional holes in 
the simplicial complex. 

Applications of alpha shapes have as yet focussed on biological systems. Their 
use in characterizing th e topology and structure of macromo l ecules. The work b y 
Liang and collaborators (lEdelsbrunner. Facello & Liang . 1 998 : Liang et al. . 1 998allbl : 



Liang. Edelsbrunner & Woodwardl 1998cl) uses alpha shapes and betti numbers to 
assess the voids and pockets in an effort to classify complex protein structures, a 
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Fig. 49. Examples of alpha shapes of the LCDM GIF simulation. Shown are central slices 
through two alpha shapes (top: low alpha; bottom: high alpha). The image shows the sen- 
sitivity of alpha shapes to the topology of the matter distribution. From: Vegter et al. 2007. 
Courtesy: Bob Eldering. 
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highly challenging task given the 10,000-30,000 protein families involving 1,000- 
4,000 complic ated folds. Given the interest in the topology of the cosmic mass dis- 
tribution (e.g. iGott. Dickinson & Melotll Il986t iMecke. Buchert & Wagneii Il994 : 
Schmalzing et al. , 19991) . it is evident that alpha shapes also provide a highly inter- 



esting tool for studying the topology of the galaxy distribution and N-body simula- 
tions of cosmic structure formation. Directly connected to the topology of the point 
distribution itself it would discard the n eed of user-defined filter kernels. 

In a recent study lVegter et al.l (l2007h computed the alpha shapes for a set of GIF 
simulations of cosmic structure formation (see fig.|49]l. On the basis of a calibration 
of the inferred Minkowski functionals and Betti numbers from a range of Voronoi 
clustering models their study intends to refine the knowledge of the topological 
structure of the Cosmic Web. 



15.6 the Multiscale Morphology Fiher 



The multiscale detection technique - MMF - is used for characterizing the differ- 
ent morphological elemen ts of the large scale matter distribution in the Universe 
Aragon-Calvo et alj ( 120071) . The method is ideally suited for extracting catalogues 
of clusters, walls, filaments and voids from samples of galaxies in redshift surveys 
or particles in cosmological N-body simulations. 

The multiscale filter technology is particularly oriented towards recognizing and 
identifying features in a density or intensity field on the basis of an assessment of 
their coherence along a range of spatial scales and with the virtue of providing a 
generic framework for characterizing the local morphology of the density field and 
enabling the selection of those morphological features which the analysis at hand 
seeks to study. The Multiscale Morphology Filter (MMF) method has been devel- 
oped on the basis of visuali zation and feature ex traction techniques in computer 



vision and medical research (iFlorack et al.L 119921) . The technology, finding its ori 



gin in computer vision researc h, has been optimiz ed within the context of feature 



detections in medical imaging. iFrangi et al.l (119981) and I Sato et al.l (119981) presented 



its operation for the specific situation of detecting the web of blood vessels in a med- 
ical image. This defines a notoriously complex pattern of elongated tenuous features 
whose branching make it closely resemble a fractal network. 

The density or intensity field of the sample is smoothed over a range of multiple 
scales. Subsequently, this signal is processed through a morphology response filter 
Its form is dictated by the particular morphological feature it seeks to extract, and 
depends on the local shape and spatial coherence of the intensity field. The morphol- 
ogy signal at each location is then defined to be the one with the maximum response 
across the full range of smoothing scales. The MMF translates, extends and opti- 
mizes this technology towards the recognition of the major characteristic structural 
elements in the Megaparsec matter distribution. It yields a unique framework for the 
combined identification of dense, compact bloblike clusters, of the salient and mod- 
erately dense elongated filaments and of tenuous planar walls. Figure|50]includes a 
few of the stages involved in the MMF procedure. 



Geometric Analysis Cosmic Web 



113 



Crucial for the ability of the method to identify anisotropic features such as 
filaments and walls is the use of a morphologically unbiased and optimized con- 
tinuous density field retaining all features visible in a discrete galaxy or particle 
distribution. Accordingly, DTFE is used to process the particle distribution into a 
continuous density field (top centre frame, fig.lSOt. The morphological intentions of 
the MMF method renders DTFE a key element for translating the particle or galaxy 
distribution into a representative continuous density field /dtfe- 



Structure 


A ratios 


A constraints 


Blob 
Line 
Sheet 


A[ - A2 - At, 
Al — A2 ^ At, 

Ai^ A2- As 


^3 < ; A2<0; Ai<0 
A3<0; A2<0 
^3 < 



Table 2. Behaviour of the eigenvalues for the characteristic morphologies. The lambda con- 
ditions describe objects with intensity higher that their background (as clusters, filaments and 
walls). For voids we must reverse the sign of the eigenvalues. From the constraints imposed 
by the A conditions we can describe the blob morphology as a subset of the line which is itself 
a subset of the wall. 



In the implementation of lAragon-Calvo et al.l ( 120071) the scaled representations 



of the data are obtained by repeatedly smoothing the DTFE reconstructed density 
field /dtfe with a hierarchy of spherically symmetric Gaussian filters Wq having 
different widths R: 



fs{x) = J dy /dtfe(3') Wciy, x) 

where Wc denotes a Gaussian filter of width R. A pass of the Gaussian smoothing 
filter attenuates structure on scales smaller than the filter width. The Scale Space 
itself is constructed by stacking these variously smoothed data sets, yielding the 
family of smoothed density maps /„ : 



levels 

In essence the Scale Space structure of the field is the {D + I) dimensional space 
defined by the D-dimensional density or intensity field smoothed over a continuum 
of filter scales Rq- As a result a data point can be viewed at any of the scales where 
scaled data has been generated. 

The crux of the concept is that the neighbourhood of a given point will look dif- 
ferent at each scale. While there are potentially many ways of making a comparison 
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Fig. 50. Schematic overview of the Multiscale Morphology Filter (MMF) to isolate and ex- 
tract elongated filaments (dark grey), sheetlike walls (light grey) and clusters (black dots) in 
the weblike pattern of a cosmological N-body simulation. The first stage is the translation of 
a discrete particle distribution (top lefthand frame) into a DTFE density field (top centre). The 
DTFE field is filtered over a range of scales (top righthand stack of filtered fields). By means 
of morphology filter operations defined on the basis of the Hessian of the filtered density 
fields the MMF successively selects the regions which have a bloblike (cluster) morphology, 
a filamentary morphology and a planar morphology, at the scale at which the morphological 
signal is optimal. This produces a feature map (bottom lefthand). By means of a percola- 
tion criterion the physically significant filaments are selected (bottom centre). Following a 
sequence of blob, filament and wall filtering finally produces a map of the different morpho- 
logical features in the particle distribution (bottom lefthand). The 3-D isodensity contours in 
the bottom lefthand frame depict the most pronounced features. From Aragon-Calvo et al. 
2007. 



of the scale dependence of local environment. I Aragon-Calvo et al.l (120071) chose to 
calculate the Hessian matrix in each of the smoothed replicas of the original data 
to describe the the local "shape" of the density field in the neighbourhood of that 
point. In terms of the Hessian, the local variations around a point xo of the density 
field fix) may be written as the Taylor expansion 



fixQ + s) = fixo) + «^V/(ji:o) + y'^'J-(ixo)s + ... 



(58) 



where 

fxx f\x fzx 

'H = fxy fyy f,y (59) 
^ fxz fyz fzz y 

is the Hessian matrix. Subscripts denote the partial derivatives of / with respect to 
the named variable. There are many possible algorithms for evaluating these deriva- 
tives. In practice, the Scale Space procedure evaluates the Hessian directly for a 
discrete set of filter scales by smoothing the DTFE density field by means of Mexi- 
can Hat filter, 



dxjdxj 



fsix) - /dtfe 



-I 



dxjdxj 

Ayfiy) 



(xj - ydixj - yj) - dijRl 
Ri 



WG(y,x) 



(60) 



with xi,X2,X2 = x,y,z and Sjj the Kronecker delta. In other words, the scale space 
representation of the Hessian matrix for each level « is evaluated by means of a 
convolution with the second derivatives of the Gaussian filter, also known as the 
Marr or "Mexican Hat" Wavelet. 

The eigenvalues of the Hessian matrix at a point encapsulate the information on 
the local shape of the field. Eigenvalues are denoted as Aa{x) and arranged so that 
Ai>A2> A3: 
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- Aa{x) 6ij 



= 0, fl = 1,2,3 (61) 



dxidxj 

with Ai > A2> /I3 

The Aiix) are coordinate independent descriptors of the behaviour of the density 
field in the locality of the point x and can be combined to create a variety of morpho- 
logical indicators. They quantify the rate of change of the field gradient in various 
directions about each point. A small eigenvalue indicates a low rate of change of the 
field values in the corresponding eigen-direction, and vice versa. The corresponding 
eigenvectors show the local orientation of the morphology characteristics. 

Evaluating the eigenvalues and eigenvectors for the renormalised Hessian '7Y of 
each dataset in a Scale Space shows how the local morphology changes with scale. 
First regions in scale space are selected according to the appropriate morphology 
filter, identifying bloblike, filamentary and wall-like features at a range of scales. 
The selections are made according to the eigenvalue criteria listed in table 115.61 
Subsequently, a sophisticated machinery of filters and masks is applied to assure 
the suppression of noise and the identification of the locally relevant scales for the 
various morphologies. 

Finally, for the cosmological or astrophysical purpose at hand the identified spa- 
tial patches are tested by means of an erosion threshold criterion. Identified blobs 
should have a critical overdensity corresponding to virialization, while identified 
filaments should fulfil a percolation requirement (bottom central frame). By succes- 
sive repetition for the identification of blobs, filaments and sheets - each with their 
own morphology filter - the MMF has dissected the cosmological density field into 
the corresponding features. The box in the bottom lefthand frame shows a segment 
from a large cosmological N-body simulation: filaments are coloured dark grey, the 
walls light grey and the clusters are indicated by the black blobs. 

Once these features have been marked and identified by MMF, a large vari- 
ety of issues may be adressed. An important issue is that of environmental influ- 
ences on the formation of galaxies. The MMF identification of filaments allowed 



dAragon-Calvo et al.L l2007h to show that galaxies in filaments and walls do indeed 



have a mass-dependent alignment. 



15.7 the Watershed Void Finder 



The Watershed Void Finder (WVF) is an implementation of the Watershed Trans 
form for segmentation of images of the galaxy and matter distrib ution in t o dis- 
tinct regions and objects and the subsequent identification of voids ( PlatenL 2005 



Platen, van de Weygaert & JonesL 120071) . The watershed transform is a concept de- 
fined within the context of mathematical morphology. The basic idea behind the 
watershed transform finds its origin in geophysics. It delineates the boundaries of 
the separate domains, the basins, into which yields of e.g. rainfall will collect. The 
analogy with the cosmological context is straightforward: voids are to be identi- 
fied with the basins, while iht filaments and walls of the cosmic web are the ridges 
separating the voids from each other 
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Fig. 51. Three frames illustrating the principle of the watershed transform. The lefthand frame 
shows the surface to be segmented. Starting from the local minima the surrounding basins of 
the surface start to flood as the water level continues to rise (dotted plane initially below the 
surface). Where two basins meet up near a ridge of the density surface, a "dam" is erected 
(central frame). Ultimately, the entire surface is flooded, leaving a network of dams defines 
a segmented volume and delineates the corresponding cosmic web (righthand frame). From: 
Platen, van de Weygaert & Jones 2007. 



The identification of voids in the cosmic matter distribution is hampered by 
the absence of a clearly defined criterion of what a void is. Unlike overdense and 
virialized clumps of matter voids are not genuinely defined physical objects. The 
boundary and identity of voids is therefore mostly a m atter of definition. As a con- 
sequence there is a variety of yoidfinding algorithm s ( Kauffmann & Fairalll 1991 



El-Ad. Piran & daCostalll996t|Plionis & Basilakosll2002l:lHovle & Vogele\l l200^ 



Shandarin et al 



2006; 'Col berg et aU l2005bl: IPatiri et al.L |2006|) . A recent study 



dColberg. Pearce et al.i, i2007[) contains a balanced comparison of the performance 



of the vari ous void finders with r espect to a small region taken from the Millennium 



simulation lSpringel et al.l ( 120051) 



Watersheds 

With respect to the other void finders the watershed algorithm seeks to define a 
natural formalism for probing the hierarchical nature of the void distribution in maps 
of the galaxy distribution and in N-body sir nulations of cosmic struct ure formation. 
The WVF has several advantages (see e.g. iMever & Beuchei , Il990l) . Because it is 



identifies a void segment on the basis of the crests in a density field surrounding a 
density minimum it is able to trace the void boundary even though it has a distored 
and twisted shape. Also, because the contours around well chosen minima are by 
definition closed the transform is not sensitive to local protrusions between two 
adjacent voids. The main advantage of the WVF is that for an ideally smoothed 
density field it is able to find voids in an entirely parameter free fashion. 

The word watershed finds its origin in the analogy of the procedure with that of a 
landscape being flooded by a rising level of water. Figure |5T]illustrates the concept. 
Suppose we have a surface in the shape of a landscape. The surface is pierced at 



118 Rien van de Weygaert & Willem Schaap 




Fig. 52. A visualization of several intermediate steps of the Watershed Void Finding (WVF) 
method. The top lefthand frame shows the particles of a slice in the LCDM GIF simulation. 
The corresponding DTFE density field is shown in the top righthand frame. The next, bottom 
lefthand, frame shows the resulting n-th order median-filtered image. Finally, the bottom 
righthand frame shows the resulting WVF segmentation computed on the basis of the median 
filtered image. From: Platen, van de Weygaert & Jones 2007. 



the location of each of the minima. As the waterlevel rises a growing fraction of the 
landscape will be flooded by the water in the expanding basins. Ultimately basins 
will meet at the ridges corresponding to saddlepoints in th e density field. These 
define the boundaries of the basins, enforced by means of a sufficiently high dam. 
The final result of the completely immersed landscape is a division of the landscape 
into individual cells, separated by the ridge dams. 
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Formalism 

The WVF consists of eightfold crucial steps. The two essential first steps relate di- 
rectly to DTFE. The use of DTFE is essential to infer a continuous density field from 
a given N-body particle distribution of galaxy redshift survey. For the success of the 
WVF it is of utmost importance that the density field retains its morphological char- 
acter, ie. the hierarchical nature, the weblike morphology dominated by filaments 
and walls, and the presence of voids. In particular in and around low-density void 
regions the raw density field is characterized by a considerable level of noise. In 
an essential second step noise gets suppressed by an adaptive smoothing algorithm 
which in a consecutive sequence of repetitive steps determines the median of den- 
sities within the contiguous Voronoi cell surrounding a point. The determination of 
the median density of the natural neighbours turns out to define a stable and asymp- 
totically converging smooth density field fit for a proper watershed segmentation. 

Figure|52]is an illustration of four essential stages in the WVF procedure. Start- 
ing from a discrete point distribution (top left), the continuous density field is de- 
termined by the DTFE procedure (top right). Following the natural smoothing by 
Natural Neighbour Median filtering (bottom left), the watershed formalism is able 
to identify the void segments in the density field (bottom right). 

16 Concluding Remarks 

This extensive presentation of tessellation-based machinery for the analysis of web- 
like patterns in the spatial matter and galaxy distribution intends to provide the inter- 
ested reader with a framework to exploit and extend the large potential of Voronoi 
and Delaunay Tessellations. Even though conceptually not too complex, they do 
represent an intriguing world by themselves. Their cellular geometry paves the way 
towards a natural analysis of the intricate fabric known as Cosmic Web. 

17 Acknowledgments 

RvdW wishes to thank Vicent Martinez, Enn Saar and Maria Pons for their invi- 
tation and hospitality during this fine and inspiring September week in Valencia, 
and their almost infinite patience regarding my ever shifting deadlines. We owe 
thanks to O. Abramov, S. Borgani, M. Colless, A. Fairall, T. Jarrett, M. JuriC, R. 
Landsberg, O. Lopez-Cruz, A. McEwen, M. Ouchi, J. Ritzerveld, M. Sambridge, 
V. Springel and N. Sukumar for their figures used in this manuscript. The work 
described in this lecture is the result of many years of work, involving numerous 
collaborations. In particular we wish to thank Vincent Icke, Francis Bernardeau, Di- 
etrich Stoyan, Sungnok Chiu, Jacco Dankers, Inti Pelupessy, Emilio Romano-Diaz, 
Jelle Ritzerveld, Miguel Aragon-Calvo, Erwin Platen, Sergei Shandarin, Gert Veg- 
ter, Niko Kruithof and Bob Eldering for invaluable contributions and discussions, 
covering nearly all issues touched upon in this study. In particularly grateful we are 



120 



Rien van de Weygaert & Willem Schaap 



to Bernard Jones, for his enthusiastic and crucial support and inspiration, and the 
many original ideas, already over two decades, for all that involved tessellations, 
DTFE, the Universe, and much more .... RvdW is especially grateful to Manolis and 
Menia for their overwhelming hospitality during the 2005 Greek Easter weeks in 
which a substantial fraction of this manuscript was conceived ... What better set- 
ting to wish for than the view from the Pentelic mountain overlooking the cradle 
of western civilization, the city of Pallas Athena. Finally, this work could not have 
been finished without the patience and support of Marlies ... 

References 

Abazajian K., et al. (The SDSS collaboration): Asti'on. J. 126, 2081 (2003) 
Abramov O., McEwen A.S.: Int. J. Rem. Sens. 25, 669 (2004) 
Aikio J., Mahonen P: Astrophys. J. 497, 534 (1998) 

AUiez P., Cohen-Steiner D., Yvinec M., Desbrun M.: ACM Trans. Graphics 24, 517 
(2005) 

Alliez P., Ucelli G., Gotsman M., Attene M.; Recent Advances in Remeshing of 

Surfaces. Report of the AIM@SHAPE EU network (2005) 
Amenta N., Bern M., KamvysseUs M.: Siggraph '98, 415 (1998) 
Amenta N., Bern M.: Discrete and Comp. Geometry 22, 481 (1999) 
Arad I., Dekel A., Klypin A.: Mon. Not. R. Astron. Soc. 353, 15 (2004) 
Aragon-Calvo M.A.: Morphology and Dynamics of the Cosmic Web, Ph.D. thesis, 

Groningen University (2007) 
Aragon-Calvo M.A., Jones B.J.T., van de Weygaert R., van der Hulst J.M.; Mon. 

Not. R. Astron. Soc, subm. (2007) 
Arbabi-Bidgoli S., Miiller V.: Mon. Not. R. Astron. Soc. 332, 20 (2002) 
Arnol'd V.I.: The Theory of Singularities and Its Applications, (Cambridge Univ. 

Press, Cambridge 1991) 
Ascasibar Y., Binney J.: Mon. Not. R. Astron. Soc. 356, 872 (2005) 
Babul A., Starkman G.D.: Astrophys. J. 401, 28 (1992) 
Baccelli E, Zuyev S.: Operations Research 47, 619 (1999) 
Balbus S.A., Hawley J.E: Revs. Mod. Phys. 70, 1 (1998) 
Bertschinger E.: Astrophys. J.323, L103 (1987) 

Barrow J.D., Bhavsai- S.P, Sonoda D.H.: Mon. Not. R. Astron. Soc. 216, 17 (1985) 
Basilakos S., PHonis M., Rowan-Robinson M.: Mon. Not. R. Astron. Soc. 323, 47 
(2001) 

Bennett C.L., et al.: Astrophys. J. Suppl. 148, 1 (2003) 
Bentley J.L.: Comm. ACM 18, 509 (1975) 

Bentley J.L., Friedmann J.H.: IEEE Trans. Comput. C-27, 97 (1978) 

Bernardeau, F: Astrophys. J. 390, L61 (1992) 

Bernardeau, F: Astrophys. J. 433, 1 (1994) 

Bernardeau F, Kofman L.: Astrophys. J. 443, 479 (1995) 

Bernardeau F, Juszkiewicz R., Dekel A., Bouchet F. R.: Mon. Not. R. Astron. Soc. 
274, 20(1995) 



Geometric Analysis Cosmic Web 



121 



Bernardeau E, van de Weygaert R.: Mon. Not. R. Astron. Soc. 279, 693 (1996) 
Bernardeau E, van de Weygaert R., Hivon E., PBouchet E: Mon. Not. R. Astron. 

Soc. 290, 566(1997) 
Bernardeau E, Colombi S., Gaztanaga E., Scoccimarro R.: Physics Reports 367, 1 

(2002) 

Beucher S., Lantuejoul C. In: Pmc. Intern. Workshop on Image Processing, 
(CCETT/IRISA, Rennes, Erance 1979) 

Beucher S., Meyer E: The Morphological Approach to Segmentation: The Water- 
shed Transformation. In: Mathematical Morphology in Image Processing, ed. by 
E. Dougherty (M. Dekker Inc., New York 1993) 

Beyer T., Schaller G., Deutsch A., Meyer-Hermann M.: Comput. Phys. Commun. 
172, 86 (2005) 

Beyer T, Meyer-Hermann M.: WSEAS Trans. Syst. 5, 1 100 (2006) 
Boissonnat J.-D., Yvinec M., Bronniman H.: Algorithmic Geometry (Cambr. Univ. 
Press 1998) 

Boissonnat J.-D., Cohen-Steiner D., Mourain B., Rote G., Vegter G.: Meshing of 
surfaces. In: Effective Computational Geometry for Curves and Surfaces, ed. by 
J.-D. Boissonnat, M. Teillaud (Springer- Verlag 2006) 

Boots B.N.: Metallography 17, 411 (1984) 

Borgani S., Guzzo L.: Nature 409, 39 (2001) 

Bradac M., Schneider P., Lombard! M., Steinmetz M., Koopmans L.V.E., Navarro 

J.E: Astron. Astrophys. 423, 797 (2004) 
Branchini E., Teodoro L., Erenk C.S., Schmoldt I., Efstathiou G., White S.D.M., 

Saunders W., Sutherland W., Rowan-Robinson M., Keeble O., Tadros H., Maddox 

S., Oliver S.: Mon. Not. R. Astron. Soc. 308, 1 (1999) 
Braun J., Sambridge M.: Nature 376, 655 (1995) 

Brown G.S.: New Zealand Eorestry Service Research Notes 38, 1 (1965) 
Buhmann M.D.: Approximation and interpolation with radial functions. In: Multi- 
variate Approximation and Applications, ed. by N. Dyn, D. Leviatan, D. Levin, 
A. Pinkus (Cambr Univ. Press, Cambridge 2001), pp. 25-43 
Cappellari M., Copin Y.: Mon. Not. R. Astron. Soc. 342, 345 (2003) 
Cheng S.-W., Dey T.K., Ramos E., Ray T.: Sampling and meshing a surface with 
guaranteed topology and geometry. In: Proc. 20th Sympos. Comput. Geom. 2004 
Chiang L.-Y, Coles P: Mon. Not. R. Astron. Soc. 311, 809 (2000) 
Chincarini G., Rood H.J.: Nature 257, 294 (1975) 

Chiu S.N., van de Weygaert R., Stoyan D.: Journ. Appl. Prob. 28, 356 (1996) 
Christ N.H., Eriedberg R., Lee TD.: Nuclear Physics B202, 89 (1982a) 
Christ N.H., Eriedberg R., Lee TD.: Nuclear Physics B210, 310 (1982b) 
Christ N.H., Eriedberg R., Lee TD.: Nuclear Physics B210, 337 (1982c) 
Colberg J.M., Krughoff K.S., Connolly A.J.: Mon. Not. R. Astron. Soc. 359, 272 
(2005) 

Colberg J.M., Sheth R.K., Diaferio A., Gao L., Yoshida N.: Mon. Not. R. Astron. 

Soc. 360, 216 (2005b) 
Colberg J.M.: Mon. Not. R. Astron. Soc. 375, 337 (2007) 



122 



Rien van de Weygaert & Willem Schaap 



Colberg J.M., Pearce R, Brunino R., Foster C, Platen E., Basilakos S., Fairall A., 
Feldman H., Gottlober S., Hahn O., Hoyle F, Miiller V., Nelson L., Neyrinck M., 
Plionis M., Porciani C, Shandarin S., Vogeley M., van de Weygaert R.: Mon. Not. 
R. Astron. Soc, subm. (2007) 

Coles P, Chiang L.-Y.: Nature 406, 376 (2000) 

Colless M., et al.: aitto-ph/0306581 (2003) 

Colombi S., Pogosyan D., Souradeep T.: Phys. Rev. Let. 85, 5515 (2000) 
Cressie N.: Statistics for Spatial Data, rev. ed. (John Wiley & Sons, Chichester 
1993) 

de Berg M., van Kreveld M., Overmars M., Schwarzkopf O.: Computational Geom- 
etry: Algorithms and Applications, 2nd ed. (Springer- Verlag Heidelberg 2000) 

Delaunay B. N.: Bull. Acad. Sci. USSR: Clase Sci. Mat. 7, 793 (1934) 

Dey T., Edelsbrunner H., Guha S.: Computational Topology. In Adances in Dis- 
crete and Computational Geometry, ed. by B. Chazelle, J.E. Goodman, R. Pollack 
(Contemporary Mathematics, AMS, Providence 1998) 

Dey T.K.: Curve and surface reconstruction. In: Handbook of Discrete and Com- 
putational Geometry, 2nd. ed., ed. J.E. Goodman, J. O'Rourke (CRC Press LLC, 
Baton Raton 2004) 

Dinh H.Q., Turk G., Slabaugh G.: IEEE Trans. Pattern Anal. Mach. Intell. 24, 1358 
(2002) 

Dirichlet G.L.: Journal fur die Reine und Angwandte Mathematik 40, 209 (1850) 
Dubinski J., da Costa L.N., Goldwirth D.S., Lecar M., Piran T.: Astrophys. J. 410, 
458(1993) 

Ebeling H., Wiedenmann G.: Phys. Rev. E 47, 704 (1993) 

Edelsbrunner H., Kirkpatrick D., Seidel R.: IEEE Trans. Inform. Theory IT-29, 55 1 
(193) 

Edelsbrunner H., Miicke E.P: ACM Trans. Graphics 13, 43 
Edelsbrunner H., ShahN.R.: Algorithmica 15, 223 (1996) 
Edelsbrunner H., Facello M., Liang J.: Discrete Appl. Math. 88, 83 (1998) 
Edelsbrunner H., Letscher D., Zomorodian A.: Discrete and Computational Geom- 
etry 28, 511 

Einasto J., Joeveer M., Saar E.: Nature, 283, 47 (1980) 

El-Ad H., Piran T., da Costa L.N.: Astrophys. J. 462, L13 (1996) 

El-Ad H., Piran T.: Astrophys. J. 491, 421 (1997) 

Eldering B.: Topology of galaxy models. M.Sc. thesis, Groningen University (2006) 
Espanol P, Phys. Rev. E. 57, 2930 (1998) 

Euler L.: Letter dated 1759 October 27 from Euler to Lagrange, Opera Postuma 2, 

559(1862) 

De Fabritiis G., Coveney PV.:'arXiv:cond-mat/ 0301378| (2003) 

Fairall A. P.: Voids in the 6dF Galaxy Survey. In: Cosmic Voids: Much ado about 

Nothing, ed. R. van de Weygaert, M. Vogeley, R. Sheth, P.J.E. Peebles, F. Hoyle 

(Edita 2008) 

Fischer R.A., Miles R.E.: Mathematical Biosciences 18, 335 (1973) 



Geometric Analysis Cosmic Web 



123 



Fisher K. B., Huchra J. P., Strauss M. A., Davis M., Yahil A., Schlegel D.: Astro- 
phys. J. Suppl. 100, 69(1995) 

Flekk0y E.G., Coveney P.V., De Fabritiis G.: Phys. Rev. E. 62, 2140 (2000) 

FlorackL.M.J., ter Haar Romeny B.M., Koenderink J.J., ViergeverM.A.: Image and 
Vision Computing 10, 376 (1992) 

Fortune S.: Voronoi Diagrams and Delaunay Trinagulations. In: Computing in Eu- 
clidian Geometry, ed. D.-Z. Du, F. Hwang (World Scientific, Lect. Notes Ser. 
Computing 1, 1992) 

Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A.; Multiscale Vessel En- 
hancement Filtering. In: Medical Image Computing and Computer-Assisted In- 
terventation - MICCAI'QS, Lecture Notes in Computer Science 1496 (Springer- 
Verlag, Beriin Heidelberg 1998) 

Friedman J.H., Bentley J.L., Finkel R.A.: ACM Trans. Math. Softw. 3, 209 (1977) 

Fukugita M., Peebles PJ.E.: Astrophys. J. 616, 643 (2004) 

Gaztanaga E.: Astrophys. J. 398, L17 (1992) 

Geller M., Huchra J. P: Nature 246, 897 (1989) 

Gerald C. F, Wheatley P.O.: Applied Numerical Analysis, (Addison Wesley, Read- 
ing 1999) 

Gilbert E.N.: Annals Math. Statistics 33, 958 (1962) 

Gingold R.A., Monaghan J.J.: Mon. Not. R. Astron. Soc. 181, 375 (1977) 

Gnedin N.Y.: Astrophys. J. Suppl. 97, 231 (1995) 

Gold CM., Remmele PR., Roos T.: Voronoi methods in GIS. In: Lecture Notes in 
Computer Science 1340 (Springer- Verlag Berlin 1997) pp. 21-35 

Goodman J. E., O'Rourke J. (editors): Handbook of Discrete and Computational 
Geometry, 2nd ed. (Chapman & Hall/CRC press 1980) 

Gor v., Shapiro B.E., Jonsson H., Heisler M., Reddy G. V, Meyerowitz E.M., Mjol- 
sness E.: A Software Architecture for Developmental Modelling in Plants: The 
Computable Plant Project. In: Bioinformatics of Genome Regulation and Struc- 
ture n (Springer US 2006) 

Gott J.R. Ill, Dickinson M., Melott A.L.: Astrophys. J. 306, 341 (1986) 

Gott J.R. Ill, Juric M., Schlegel D., Hoyle F, Vogeley M., TegmarkM., BahcallN., 
Brinkmann J.: Astrophys. J. 624, 463 (2005) 

Gregory S.A., Thompson L.A.: Astrophys. J. 222, 784 (1978) 

Hardy R.L.: Journal of Geophysical Research 76, 1905 

Hayashino T. et al.: Asti'on. J. 128, 2073 (2004) 

Heath Jones D. et al.: Mon. Not. R. Astron. Soc. 355, 747 (2004) 

Hoffman Y., Ribak E.: Astrophys. J. 380, L5 (1991) 

Hoyle F, Vogeley M.: Astrophys. J. 566, 641 (2002) 

Hoyle F, Vogeley M.: Astrophys. J. 580, 663 (2002) 

Icke v., van de Weygaert R.: Astron. Astrophys. 184, 16 (1987) 

IlUev I.T., Ciardi B., Alvarez M.A., MaselH A., FeiTara A., Gnedin N.Y., Mellema 
G., Nakamoto M.L., Norman M.L., Razoumov A.O., Rijkhorst E.-J., Ritzerveld 
J., Shapiro PR., Susa H., Umemura M., Whalen D.J.: Mon. Not. R. Astron. Soc. 
371, 1057 (2006) 



124 



Rien van de Weygaert & Willem Schaap 



Iske A.: Scattered Data Modelling Using Radial Basis Functions. In: Tutorials 
on Multiresolution in Geometric Modelling, ed. A. Iske, E. Quak, M.S. Floater 
(Mathematics and Visualization, Springer- Verlag, Heidelberg 2002) pp. 287-315 

Jarrett T.H.: PAS A 21, 396 (2004) 

Jones B.J.T., Martinez V.J., Saar E., Trimble V.: Rev. Mod. Phys. 76, 1211 (2005) 
Juszkiewicz R., Weinberg D.H., Amsterdamski P., Chodorowski M., Bouchet F.: 

Astrophys.J. 442, 39(1995) 
Kaiser N.: Statistics of Gravitational Lensing 2: Weak Lenses. In: New Insights into 

the Universe, Lecture Notes in Physics 408, ed. by V.J. Martinez, M. Portilla, D. 

Saez (Springer- Verlag Berlin Heidelberg New York 1992), pp. 279 
Kaiser N., Squires G.: Astrophys. J. 404, 441 (1993) 

Kansal A.R., Torquato S., Harsh G.R., Chiocca E.A., Deisboeck T.S.: J. Theor Biol. 
203, 367 

Katz N., Weinberg D.H., Hernquist L..: Astrophys. J. Suppl. 105, 19 (1996) 
Kauffmann G., Fairall A.P: Mon. Not. R. Astron. Soc. 248, 313 (1991) 
Kauffmann G., Colberg J.M., Diaferio A., White S.D.M.: Mon. Not. R. Astron. Soc. 

303, 188 (1999) 
Kendall D.G.: Stat. Science 4, 87 (1989) 
Kiang T.: Z. Astrophys. 64, 433 (1966) 

Kim R. S. J., Kepner J.V, Postman M., Strauss M.A., Bahcall N.A., Gunn J.E., 
Lupton R.H., Annis J., Nichol R.C., Castander F.J., Brinkmann J., Brunner R.J., 
Connolly A., Csabai I., Hindsley R.B., Izevic Z., Vogeley M.S., York D.G.: As- 
tron. J. 123, 20 (2002) 

Kimia B.B., Leymarie F.F.: Symmetry-based representations of 3D data. In: Proc. 
2001 Intern. Conf. on Image Processing 2, 581 (2001) 

Kirshner R.P, Oemler A., Schechter PL., Shectman S.A.: Astrophys. J. 248, L57 
(1981) 

Kirshner R.P, Oemler A., Schechter PL., Shectman S.A.: Astrophys. J. 314, 493 
(1987) 

Koumoutsakos P: Ann. Rev. Fluid Mech. 37, 457 (2005) 

Lachieze-Rey M., da Costa L.N., Maurogordata S.: Astrophys. J. 399, 10 (1992) 

Lagrange J.L.: Oeuvres de Lagrange 1, 151 (1762) 

Lancaster P., Salkauskas K.: Mathematics of Computation 37, 141 (1981) 

de Lapparent V, Geller M.J., Huchra J.P: Astrophys. J. 302, LI (1986) 

de Lapparent V, Geller M.J., Huchra J.P: Astrophys. J. 369, 273 (1991) 

Lawson C.L.: Software for C' surface interpolation. In: Mathematical Software III, 

vol. 3, ed. by J.R. Rice (Academic Press, New York 1977) 
Leymarie FF, Kimia B.B.: Computation of the shock scaffold for unorganized point 
clouds in 3D. In: Proc. 2003 Comp. Soc. Conf. on Computer Vision and Pattern 
Recognition 1, 1-821 (2003) 
Li G.-L., Mao S., Jing YP, Kang X., Bartelmann M.: Astrophys. J. 652, 43 (2006) 
Liang J., Edelsbrunner H., Fu P., Sudhakar PV, Subramaniam S.: Proteins: Struc- 
ture, Function, and Genetics 33, 1 (1998a) 



Geometric Analysis Cosmic Web 



125 



Liang J., Edelsbrunner H., Fu P., Sudhakar P.V., Subramaniam S.: Proteins: Struc- 
ture, Function, and Genetics 33, 18 (1998b) 
Liang J., Edelsbrunner H., Woodward C: Protein Science 7, 1884 (1998) 
Lodha S.K., Franke R.: Scattered Data Techniques for Surfaces. In: Pwc. Dagstuhl 

Conf. Scientific Visualization, (IEEE Computer Society Press), pp. 182-222 
Lombardi M., Schneider P: Astron. Astrophys. 373, 359 (2001) 
Lombardi M., Schneider P: Astron. Astrophys. 392, 1 153 (2002) 
Lombardi M., Schneider P: Astron. Astrophys. 407, 385 (2003) 
Lopes PA.A., de Carvalho R.R., Gal R.R., Djorgovski S.G., Odewahn S.C., Maha- 

bal A.A., BrunnerR.J.: Astron. J. 128, 1017 (2004) 
Lucy L.B.: Astron. J. 82, 1013 (1977) 
Luo R., Vishniac E.: Astrophys. J. Suppl. 96, 429 (1995) 

Maddox S.J., Efstathiou G., Sutherland W.J., Loveday J.: Mon. Not. R. Astron. Soc. 
242, 43 (1990a) 

Maddox S.J., Efstathiou G., Sutherland W.J., Loveday J.: Mon. Not. R. Astron. Soc. 
243, 692 (1990b) 

Maddox S.J., Efstathiou G., Sutherland W.J.: Mon. Not. R. Astron. Soc. 246, 433 
(1990c) 

Marinoni C, Davis M., Newman J.A., Coil A.L.: Astrophys. J. 580, 122 (2002) 
Martinez V., Saar E.: Statistics of the Galaxy Distribution, (Chapman & Hall/CRC 
press 2002) 

Martinez V., Starck J.-L., Saar E., Donoho D.L., Reynolds S.C., de la Cruz P., Pare- 

des S.: Astrophys. J. 634, 744 (2005) 
Mavripilis D.J.: Ann. Rev. Fluid Mech. 29, 473 (1997) 
Mecke K.R., Buchert T., Wagner H.: Astron. Astrophys. 288, 697 (1994) 
Meijering J.L.: Philips Research Reports 8, 270 (1953) 
Mellier Y.: Ann. Rev. Astron. Astrophys. 37, 127 (1999) 
Metropolis N., Ulam S.: J. Am. Stat. Assoc. 44, 247 (1949) 
Meyer F, Beucher S.: J. Visual Comm. Image Rep. 1, 21 (1990) 
Miles R.E.: Math. Biosciences 6, 85 (1970) 
Miles R.E.: Adv. Appl. Prob. (Suppl.) 4, 243 (1972) 

Miles R.E.: A synopsis of 'Poisson Flats in Euclidian Spaces'. In: Stochastic Geom- 
etry, ed. by E.F Harding, D.G. Kendall (John Wiley, New York 1974) pp. 202-227 
Miles R.E., Maillardet R.J.: J. Appl. Prob. 19A, 97 (1982) 
Miller W.: Class. Quant. Grav. 14, L199 (1997) 

Mitra S.K., Sicuranza G.L.: Nonlinear Image Processing, (Academic Press Inc., 
U.S. 2000) 

Molchanov S.A., Surfailis D., Woyczynski W.A.: Annals Appl. Prob. 7, 200 (1997) 
M0ller J.: Adv. Appl. Prob. 21, 37 (1989) 

M0ller J.: Lectures on Random Voronoi Tessellations. In: Lecture Notes in Statistics, 
87 (Springer- Verlag, New York 1994) 

Monaghan J.J.: Ann. Rev. Astron. Astrophys. 30, 543 (1992) 

Miicke E.P: Shapes and Implementations in three-dimensional geometry, Ph.D. the- 
sis, Univ. Illinois Urbana-Champaign (1993) 



126 



Rien van de Weygaert & Willem Schaap 



NeyrinckM. C, Gnedin N. Y, Hamilton A. J. S.: Astrophys. J. 356, 1222 (2005) 
NovikovD., Colombi S., Dore O.: Mon. Not. R. Astron. Soc. 366, 1201 (2006) 
Nusser A., Dekel A., BertschingerE., Blumenthal G. R.: Astrophys. J. 379, 6 (1991) 
Okabe A., Boots B., Sugihara K., Chiu S. N.: Spatial tessellations : concepts and ap- 
plications ofvoronoi diagrams, 2nd ed. (John Wiley &. Sons, Chichester, Toronto 
2000) 

Ord J.K.: Mathematical Scientist 3, 23 (1978) 

Patiri S.G., Betancort-Rijo J.E., Prada R, Klypin A., Gottlober S.: Mon. Not. R. 

Astron. Soc. 369, 335 (2006a) 
Patiri S.G., Prada E, Holtzman J., Klypin A., Betancort-Rijo J.E.: Mon. Not. R. 

Astron. Soc. 372, 1710 (2006b) 
Peacock J.A., Dodds S.J.: Mon. Not. R. Astron. Soc. 267, 1020 (1994) 
Peebles P. J. E.: The large-scale structure of the universe, (Princeton University 

Press 1980) 

Pelupessy E.I., Schaap W.E., van de Weygaert R.: Astron. Astrophys. 403, 389 
(2003) 

Pen Ue-li: Astrophys. J. Suppl. 115, 19 (1998) 
Pimbblet K.A.: Mon. Not. R. Astron. Soc. 358, 256 (2005) 
Platen E.: Segmenting the Universe. M.Sc. thesis, Groningen University (2005) 
Platen E., van de Weygaert R., Jones B.J.T.: Mon. Not. R. Astron. Soc, in press 
(2007) 

Plionis M., Basilakos S.: Mon. Not. R. Astron. Soc. 330, 399 (2002) 

Powell M.J.D.: The theory of radial basis function approximation in 1990. In: Ad- 
vances in Numerical Analysis II: Wavelets, Subdivision, and Radial Functions, 
ed. by W.A. Light (Oxford Univ. Press, Oxford 1992), pp. 105-210 

Press W.H., Teukolsky S.A., Vetterling W.T, Flannery B.P: Numerical Recipes, 
(Cambridge Univ. Press, Cambridge U.K. 1992) 

Ramella M., Boschin W., Fadda D., Nonino M.: Astron. Astrophys. 368, 776 (2001) 

Refregier Y.: Ann. Rev. Astron. Astrophys. 41, 645 (2003) 

Regge T: Nuovo Cimentto A 19, 558 

Regos E., Geller M.J.: Asti'ophys. J. 373, 14 (1991) 

Reiprich T.H., Bohringer H.: Astron. Nachr 320, 296 (1999) 

Ripley B.D.: Spatial Statistics, (Wiley & Sons, Chichester 1981) 

Ritzerveld J.: the simplicity of transport, Ph.D. thesis, Leiden University (2007) 

Ritzerveld J., Icke V.; Phys. Rev. E. 74, 026704 (2006) 

Ritzerveld J., Pawlik A., Schaye J.: in prep. (2007) 

Rojas R.R., Vogeley M.S., Hoyle R, Brinkmann J.: Astrophys. J. 624, 571 (2005) 
Romano-Diaz E.: Probing Cosmic Velocity Flows in the Local Universe. Ph.D. the- 
sis, Groningen University (2007) 
Romano-Diaz E., van de Weygaert R.: Mon. Not. R. Astron. Soc, in press (2007) 
Rybicki G.B., Press W.H.: Astrophys. J. 398, 169 (1992) 
Ryden B.S., Gramann M.: Astrophys. J. 383, 33 (1991) 
Sambridge M., Braun J., McQueen H.: Geophys. J. Int. 122, 837 (1995) 
Sahni V., Sathyaprakash B.S., Shandarin S.F: Astrophys. J. 495, 5 (1998) 



Geometric Analysis Cosmic Web 



127 



Sato Y., Nakajima D., Atsumi H., Koller T., Gerig G., Yoshida S., Kikinis R.: IEEE 
Medical Image Analysis 2, 143 (1998) 

Saunders W., Sutherland W.J., Maddox S.J., Keeble O., Oliver S.J., Rowan- 
Robinson M., McMahon R.G., Efstathiou G.R, Tadros H., White S.D.M., Frenk 
C.S., Carraminana A., Hawkins M.R.S.: Mon. Not. R. Astron. Soc. 317, 55 (2000) 

Schaap W. E.: the Delaunay Tessellation Field Estimator Ph.D. thesis, Groningen 
University (2007) 

Schaap W. E., van de Weygaert R.: Astron. Astrophys. 363, L29 (2000) 
Schaap W. E., van de Weygaert R.: Astron. Astrophys., subm. (2007a) 
Schaap W. E., van de Weygaert R.: Astron. Astrophys., subm. (2007b) 
Schaap W. E., van de Weygaert R.: Astron. Astrophys., subm. (2007c) 
Schaap W. E., van de Weygaert R.: Astron. Astrophys., subm. (2007d) 
Schaback R., Wendland H.: Characterization and construction of radial basis func- 
tions. In: Multivariate Approximation and Applications, ed. by N. Dyn, D. Levi- 
atan, D. Levin, A. Pinkus (Cambr. Univ. Press, Cambridge 2001), pp. 1-24 
Schaller G., Meyer-Hermann M.: Comput. Phys. Commun. 162, 9 (2004) 
Schaller G., Meyer-Hermann M.: Phys Rev. E. 71, 051910.1 (2005) 
Schmalzing J., Buchert T., Melott A.L., Sahni V., Sathyaprakash B.S., Shandarin 

S.E: Astrophys. J. 526, 568 
Schmidt J.D., Ryden B.S., Melott A.L.: Astrophys. J. 546, 609 
Schoenberg I.J.: Quart. Appl. Math. 4, 45 (1946a) 
Schoenbergl.J.: Quart. Appl. Math. 4, 112 (1946b) 
Serrano M., Espanol P, Phys. Rev. E. 64, 0461 15 (2001) 

Serrano M., De Fabritiis G., Espanol P., Flekk0y E.G., Coveney P.V.: 

J. Phys. A. Math. Gen. 35, 1605 (2002) 
Shandarin S.E, Sheth J.V., Sahni V.: Mon. Not. R. Astron. Soc. 353, 162 (2004) 
Shandarin S., Feldman H.A., Heitmann K., Habib S.: Mon. Not. R. Astron. Soc. 

376, 1629 (2006) 

Sharma S., Steinmetz M.: Mon. Not. R. Astron. Soc. 373, 1293 (2006) 

Shepard D.: A two-dimensional interpolation function for irregularly spaced points. 

In: ACM National Conference pp. 517-524 (1968) 
Sheth J. v., Sahni V., Sathyaprakash B.S., Shandarin S.E: Mon. Not. R. Astron. Soc. 

343, 22 (2003) 

Sheth R. K., van de Weygaert R.: Mon. Not. R. Astron. Soc. 350, 517 (2004) 
Shewchuk J.R.: Comp. Geom. Theor. Appl. 22, 21 (2002) 
Sibson R.: Math. Proc. Cambridge Phil. Soc. 87, 151 (1980) 

Sibson R.: A brief description of natural neighbor interpolation. In: Interpreting 

Multi-variate Data, ed. by V. Barnett (Wiley, Chichester 1981) pp. 21-36 
Sloan S.W.: Adv. Eng. Software 9, 34 
Smoot G.R, et al.: Astrophys. J. 396, LI (1992) 
Soneira R.M., Peebles PJ.E.: Astrophys. J. 211, 1 (1977) 

Sousbie T., Pichon C, Courtois H., Colombi S., Novikov D.: astroph/0602628 

(2006) 

Spergel D.N., et al. : |astro-ph/0603449K 2006) 



128 



Rien van de Weygaert & Willem Schaap 



Springel v., Yoshida N., White S.D.M.: New Astron. 6, 79 (2001) 
Springel V.: Mon. Not. R. Astron. Soc. 364, 1105 (2005) 

Springel V., White S.D.M., Jenkins A., Frenk C.S., Yoshida N., Gao L., Navarro J., 
Thacker R., Croton D., Helly J., Peacock J. A., Cole S., Thomas P., Couchman H., 
Evrard A., Colberg J.M., Pearce E: Nature 435, 629 (2005) 

Stoica R.S., Martinez V.J., Mateu J., Saar £.: Astron. Astrophys. 434, 423 (2005) 

Stoyan D., Kendall W.S., Mecke J.: Stochastic Geometry and its Applications, 1st 
ed. (John Wiley, New York 1987) 

Stoyan D., Kendall W.S., Mecke J.: Stochastic Geometry and its Applications, 2nd 
ed. (John Wiley, New York 1995) 

Sukumar N.; The Natural Element Method in Solid Mechanics. Ph.D. thesis. North- 
western University (1998) 

Szapudi 1.: Astrophys. J. 497, 16 (1998) 

TegmarkM., et al.: Astrophys. J. 606, 702 (2004) 

Thiessen A.H.: Monthly Weather Review 39, 1082 (191 1) 

Torquato S.: Random Heterogeneous Materials: Microstnicture and Macroscopic 

Properties (Springer, New York 2002) 
Turk G., O'Brien J.E: Shape transformation using variational implicit functions. 

In: Computer Graphics. Annual Conference Series (SIGGRAPH 1999) 1999, 

pp. 335-342 

Turk G., O'Brien J.F.: Reconstructing Surfaces Using Anisotropic Basis Functions. 
In: Proc. Intern. Conference Computer Vision (ICCV), Vancouver 2001, pp. 606- 
613 

Ubertini S., Succi S.: Prog. Comp. Eluid Dyn. 5, 85 (2005) 

van de Weygaert R.: The Multidimensional Binary Tree and its Use in Astronomy. 

(minor) M.Sc. Thesis, Leiden University (1987) 
van de Weygaert R.: Voids and the Geometry of Large Scale Structure. Ph.D. thesis, 

Leiden University (1991) 
van de Weygaert R.: Astron. Astrophys. 283, 361 (1994) 

van de Weygaert R.: Eroth Across the Universe, Dynamics and Stochastic Geome- 
try of the Cosmic Eoam. In: Modem Theoretical and Observational Cosmology, 
proc. 2nd Hellenic Cosmology Meeting, ed. by M. Plionis, S. Cotsakis, ASSL 
276 (Kluwer, Dordrecht 2002) pp. 1 19-272 

van de Weygaert R.: Astron. Astrophys., subm. (2007a) 

van de Weygaert R.: Astron. Astrophys., subm. (2007b) 

van de Weygaert R., van Kampen E.: Mon. Not. R. Astron. Soc. 263, 481 (1993) 
van de Weygaert R., Bertschinger E.: Mon. Not. R. Astron. Soc. 281, 84 (1996) 
Vegter G.: Computational topology. In: Handbook of Discrete and Computational 
Geometry, 2nd. ed., ed. J.E. Goodman, J. O'Rourke (CRC Press LLC, Baton Ra- 
ton 2004) pp. 719-742 
Vegter G., van de Weygaert R., Platen E., Kruithof N., Eldering B.: Mon. Not. R. 

Astron. Soc, in prep. (2007) 
Voronoi G.: J. reine angew. Math. 134, 167 (1908) 



Geometric Analysis Cosmic Web 



129 



Watson D.F.: Contouring: A Guide to the Analysis and Display of Spatial Data, 

(Pergamom Press, Oxford 1992) 
Wendland H.: Scattered Data Approximation, (Cambridge Mon. Appl. Comput. 

Math., Cambridge Univ. Press, Cambridge 2005) 
White S.D.M.: Mon. Not. R. Astron. Soc. 186, 145 (1979) 
Whitehurst R.: Mon. Not. R. Astron. Soc. 277, 655 (1995) 
Yahil A., Strauss M. A., Davis M., Huchra J. P: Astrophys. J. 372, 380 (1991) 
Zaroubi S., Hoffman Y., Fisher K.B., Lahav O.: Astrophys. J. 449, 446 (1995) 
Zel'dovich Ya. B.: Astron. Astrophys. 5, 84 (1970) 
Zel'dovich Ya.B., Einasto J., Shandarin S.F.: Nature 300, 407 (1982) 
Zomorodian A.J.: Topology for Computing, (Cambr. Mon. Appl. Comp. Math. 

2005) 



